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ABSTRACT 


(Distribution  Limitation  Statement  No.  5) 

Progress  accomplished  under  Air  Force  Contract  No.  F29601-68- 
C-0012  toward  predicting  the  electromagnetic  pulse  produced  by  a  nuclear 
weapon  is  reported  herein.  This  work  consists  of  the  following:  (1)  a  real¬ 
time  finite  difference  code,  using  prolate  spheroidal  coordinates,  which  will 
predict  the  electromagnetic  fields  produced  by  certain  current  distributions 
with  azimuthal  symmetry,  (2)  a  retarded  time  code  which  in  other  respects 
is  similar  to  the  one  listed  above,  (3)  a  one -dimensional  finite  difference 
code  for  predicting  the  fields  of  a  current  distribution  with  spherical  sym¬ 
metry,  (4)  a  computer  code  which  numerically  evaluates  the  exact  solution 
of  Maxwell's  equations  for  the  case  of  spherical  symmetry,  (5)  a  numerical 
solution  of  the  diffusion  approximation  of  Maxwell's  equations  in  a  region 
which  includes  an  infinitely  conducting  earth,  (6)  work  toward  obtaining  a 
Green's  function  for  Maxwell's  equations  with  some  restrictions  on  the  time 
dependence  of  the  conductivity,  and  (7)  a  numerical  integration  code  that 
combines  the  results  of  two  previously  existing  Monte  Carlo  neutron  and 
gamma  ray  transport  codes  in  order  to  obtain  currents  produced  by  gamma 
ray  sources  in  the  ground. 
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I.  EMP  FINITE  DIFFERENCE  CODES 


1.  INTRODUCTION 


Two  digital  computer  codes,  which  calculate  the  electromagnetic 
pulse  generated  by  a  low-altitude  nuclear  burst  in  the  presence  of  a  finitely 
conducting  earth,  have  been  developed.  It  is  assumed  in  the  codes  that  the 
currents  and  conductivities  are  known  functions  of  space  and  time,  which 
have  azimuthal  symmetry.  With  this  input.  Maxwell's  equations  are  solved 
in  prolate  spheroidal  coordinates  by  a  finite  difference  technique.  The 
geometry  for  this  coordinate  system  is  described  in  Section  1-4.  The  codes 
differ  in  that  one  is  written  in  real  time  whereas  the  other  employs  the  re¬ 
tarded  time  of  the  burst  point.  The  primary  advantages  of  the  real  time  code 
are:  (l)a  definite  physical  boundary  condition  exists  at  the  wave  front  in 
space  (i.  e.,  the  fields  are  zero  for  r>ct),  and  (2)  it  is  possible  to  have  the 
mesh  move  along  with  the  wave  front  and  thus  save  computer  storage  if  it  is 
not  necessary  to  calculate  fields  for  great  distances  behind  the  wave  front. 
Advantages  of  the  retarded  time  code  which  at  present  are  considered  de¬ 
cisive  are:  (1)  regridding  both  at  the  ground-air  interface  and  at  the  wave 
front  is  much  simpler  than  for  the  real-time  code,  and  (2)  the  spatial  deriv¬ 
atives  at  the  wave  front  in  the  air  are  smaller  than  for  the  real  time  case. 

For  these  reasons  a  production  version  of  the  retarded-time  code  is  presently 
being  produced. 

The  current  and  conductivities  which  are  now  being  used  in  the  codes 
are  calculated  from  the  equations 


J(r,  r) 


exp [a( r-r  )  -  r/R  1 

P  yJ  a 

1  +  exp  |Xor+/3  )(  r-  )J 


(1) 


-1- 


(2) 


t) 


a  exp  cAt-t  )  r/R 

_£  L  ,  P _ 1 

p2  1  +exp  J^(a'+|3)(T-Tp) 


It  is  anticipated  that  future  work  will  include  writing  codes  to  generate  more 
realistic  sources. 

Identical  problems  have  been  run  with  these  two  codes.  In  regions 
where  the  fields  have  appreciable  values  the  results  compare  quite  favorably. 
Examples  are  depicted  in  the  graphs  in  Appendix  A.  In  addition,  results  of 
these  codes  have  been  compared  to  those  of  the  one  dimensional  code  des¬ 
cribed  in  Section  1-5.  In  regions  where  the  fields  of  the  two  dimensional  codes 
are  not  affected  by  reflections  from  the  ground,  they  give  essentially  the  same 
output  as  the  one  dimensional  code  does. 

The  basic  codes  are  described  below.  Versions  are  also  written 
which  contain  regridding  at  the  ground  and  also  at  the  wave  front  for  the 
retarded  time  code. 


2.  REAL-TIME  EMP  CODE  FOR  LOW  ALTITUDE  NUCLEAR  BURSTS 
a.  Maxwell's  Equations 

In  M. K.S.  units  Maxwell's  equations  are 


V  •  D  =  p 

(3) 

V  •  B  =  0 

(4) 

Iff 

CD 

1 

(5) 

-  -  0D 

St 

(6) 

The  divergence  equations  may  be  treated  as  initial  conditions  since  Maxwell's 
equations  predict  that  if  they  are  initially  satisfied  they  will  remain  so.  Our 
attention  will  thus  be  focused  on  the  curl  equations.  In  prolate  spheroidal 
coordinates  they  are 
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We  assume  the  constitutive  equations 


D  =  eE,  B  =  pH 


(9) 


The  problem  considered  has  azimuthal  symmetry.  Thus  no  physical  quantity 
is  a  function  of  4>.  By  using  Eq.  (9)  and  this  symmetry  condition,  the  curl 
equations  are  reduced  to 
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where 
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(12) 


In  scalar  form,  Eqs.  (10)  and  (11)  separate  into  two  independent  sets  of 
coupled  equations.  One  set  involves  the  variables  E^  ,  E^  ,  ,  and 

.  The  other  set  relates  ,  B^  ,  and  j^.  Maxwell's  divergence  equa¬ 

tions  behave  similarly.  We  assume  that  for  our  problem  j ,  =  0.  It  then 

follows  that  if  E, ,  B_  ,  and  B„  are  initially  zero,  they  will  remain  so.  Thus 
9  ?  s 

we  set 


=  B?  =j<j>=° 


(13) 


The  curl  equations  become 
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In  order  to  simplify  these,  the  fields  will  be  transformed  according  to 


(16) 


(17) 


(18) 


in  which  the  primed  quantities  are  those  that  have  been  used  previously.  In 
component  form,  the  field  equations  then  become 
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Make  the  substitution 


3— »  3  +  ctE 


(22) 


Henceforth  j  will  represent  only  that  portion  of  the  current  resulting  from 
Compton  electrons. 

Finally,  we  have 
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This  set  of  equations  will  be  solved  numerically  below  by  a 
finite  difference  technique. 

b.  Difference  Equations 
Let 


5=l+(l-i)AC,  ?  =  a  +  ( j  - 1  )A?  ,  t  =  kAt  (26) 

so  that  at  a  mesh  point  a  function  of  ?  ,  ?  ,  and  t  may  be  labeled  f(  i,  j,  k ). 
Maxwell's  curl  equations  (23)  through  (25)  will  be  replaced  by  difference 
equations  centered  at  ( i,  j,  k- 1/2  ).  For  mesh  points  in  the  air,  we  find 
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Provided  the  values  of  the  fields  are  known  on  the  z  axis,  on  the  line  j-l,  and 
for  values  of  ?  greater  than  1  +  (l-i)AS  then  the  above  are  three  equations 
that  can  be  solved  simultaneously  for  the  three  fields  ( i,  j,  k  ),  ( i,  j,  k  ), 

and  B^(i,  j,  k).  From  the  symmetry  of  the  problem  it  is  easy  to  show  that  on 
the  z  axis  for  z>a, 


lim  E  { €  ,  £  ,  t )  =  0  lim  B,(?,?,t)=0  (30) 
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The  value  of  E^  at  5 
i.  e. , 


1  is  obtained  by  leapfrogging  the  field  forward  in  time. 
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Similarly,  for  £  =  a,  ?>0,  we  have 
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After  these  field  values  are  found  on  the  z  axis,  equations  (27),  (28),  and  (29) 
are  used  recursively  to  obtain  the  fields  in  the  air  mesh  on  a  line  of  constant 
J- 
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The  difference  equations  in  the  ground  are 
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The  nonzero  electric  fields  on  the  z  axis  in  the  ground  are  obtained  from 
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and  from  Eq.  (32)  by  inserting  parameters  appropriate  to  the  ground.  In  the 
above  equation,  i  is  the  value  of  i  corresponding  to  5  =  -1.  Equations  (33), 

(34),  and  (35)  are  then  employed  to  obtain  field  values  in  the  ground  along  a 
line  of  constant  i  in  a  manner  analogous  to  that  described  above  for  the  air 
mesh. 

To  treat  the  air-ground  boundary,  we  define  two  lines  of  points  labeled 

i  and  i  for  which  ?  =0.  Even  though  these  points  have  the  same  £  coordinates, 
a  S 

the  first  set  is  assumed  to  be  in  air  and  the  other  in  the  ground.  The  mesh 

points  near  the  ground  are  labeled  as  shown  in  Figure  1. 

The  appropriate  boundary  conditions  are  that  E„,  and  B,  are  con- 

?  9 

tinuous.  This  must  hold  for  all  time.  Thus  we  infer  that  3E^  /9t  and 

3B, /3 1  must  be  continuous.  Whence 
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<P  g _ 

at 


(38) 


At  [E?  (V  3’  k)_E?(la'  3’  k_1  "  8A?  [B0(lg*  3+1’  k_1  ) 

-  B.(i  +1,  j+1,  k-1  )  +  B  (i  ,  3,  k-1  )  -  B.(i  +1,  j,  k-1  ) 

9  g  9  g  9  g 

+  VV  j,  k)-B^(ig+l,  j,  k)  +  B^(ig,  j-l,  k)-B^(ig+l,  j-l,  k )] 


k<ia“l’  j+3'’  ^-1  )  -  B^i^,  j+1,  k-1  )  +  B^  ( i^-1,  j,  k-1) 


8  A?  [_9  a 


B^i^,  j,  k-1  )  +  B^(ia~l,  j,  k)  -  B^(i^,  j,  k)  +  B^d^+l,  j-l,  k) 


9  a 


JL  '  A  9 

0  a 


<T  a 


a(i  +1/2,  j,  k-1/2) 

k)]~  ~g— S7 - 


[E? ( V  ^  k-1 ) 


+  E?  (ig,  j,  k)  +  E?  (ig+l,  j>  k)  +  E?  (i  +1,  j>  k-1  )j 


a(i  -1/2,  j,  k-1/2) 
a 

— 


[E?(ia-1’  3’  k)+E?(ia'  3*  k )  +  E  ( i&,  j,  k-1) 


+  VV1'  >'  k-U]-2T  j?(ig+1/2'  S'  k-l/2)-2ic. 


j?(i&-l/2,  j,  k-1/2) 


(39) 
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2  2 

AT  [B4> ( ia’  k  ]  '  B<t> ( V  j’  k_1  }]  =  8A?  ?2  ?_  a2A?2/4 


[e  (i  +1,  j+1,  k-l)-E  (i  +1,  j,  k-l)+E  <i  +1,  j,  k) 

L  s  g  s  S  **  S 


-  E_  (i  +1,  j-1,  k)  +  E_  (i  ,  j+1,  k-1  )  -  E  ( i  ,  j,  k-1  ) 
5  g  5  g  5  g 


+  Ec  (i_,  j,  k)  -  Ec(i^  j-1,  k)  +  E?  (i&,  j+1,  k-1  )  -  Ep  (  iQ,  j,  k-1  ) 


€'VJ  «  g 


€  a 


+  E„  (i  ,  j,  k)  -  E  (i  ,  j-1,  k)  +  Ep  (i  -1,  j+1,  k-1  )  -  E  (i  -1,  j,  k-1  ) 
s  a  $  a  ?  a  s  a 


+  E  (i  -1,  j,  k)  -  E  (i  -1,  j-1,  k)l  -  -j — 4^ 

?  a  l  a  J  8A?  ^2  -a2A?2/4 


[e  „(i  -1,  j,  k)  +E  (i  -1,  j,  k-1  )  -  E  (i  +1,  j,  k-1  ) 
L  S  &  S  a  s  & 


E?(i  +1,  j,  k)+E?(ia-l,  j+1,  k-1 )  +  E?(ia-1,  j-1,  k) 


-E?(i  +1,  j+1,  k-1)  -  E?(ig+1,  3 - 1 ,  k )] 


V  g 


Evaluate  9E  /9t  at  (i  ,  j,  k-1/ 2) 

S  3- 


j,  k)  -  E  (i  ,  j,  k-l 
s  3- 


»] 


j+l ,  k-l  ) 


a(i  ,  j,  k-1/2) 

-  VV  j*  k"1)  +  B^(ia'  j'  k)  "  W  j-1'  k]-“  a 


2  €„ 


[E  (i  ,  j,  k)  +  E  (i  ,  j,  k-l)l  -  ~  j,.  (i  »  j,  k-1/2) 
L?a  £a  -I  en  €  a 


(41) 


Evaluate  3E_  / at  at  (i  ,  j,  k-1/2  ) 
?  g 


ir  [Es ( V  j'  k) "  Es ( V  j- k_1  ’] =  tk?  [V  v  j+l- k_1  ’ 

<f(i  ,  j,  k-1/2) 

-  W  j'  k~1)  +  VV  j'  k)'VV  +1’  k  !]  '  g~2'r - 


[Ef(i  ,  j,  k)+E5(ig,  j,  k-l)] j,  k-1/2) 


(42) 


After  computing  fields  in  the  air  and  ground  meshes  for  given  i  and  k,  Eqs. 
(39)  through  (42)  can  be  solved  simultaneously  to  give  the  four  independent 
field  components  at  the  boundary. 


3.  RETARDED- TIME  EMP  CODE  FOR  LOW  ALTITUDE  NUCLEAR 

BURSTS 

a.  Maxwell's  Equations 

The  time  variable  in  this  code  will  be  taken  as  the  retarded  time 
of  the  burst  point  for  waves  that  propagate  in  air.  Thus  put 


T 


(43) 
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v  ... 


t>wr  . 


where  t  is  the  retarded  time,  t  is  the  real  time,  c  is  the  velocity  of  light 
in  free  space,  and  ?  and  ?  are  described  in  Section  1-4.  The  corresponding 
transformations  for  the  derivatives  are 


(44) 


V 


a  - 

r  _a_ 
c  3t 


(45) 


Here  V  is  the  gradient  operator  at  constant  r  and  £  is  the  unit  radial  vector 
in  a  spherical  polar  coordinate  system  centered  at  the  burst  point  of  the 
weapon.  After  applying  Eq.  (43)  through  (45),  Maxwell's  equations  in  M.K.S. 
units  become 


„  -  1  a  9D 

V  •  D  -  —  r  •  — —  =  p 
1  c  9t 


(46) 


r*  1  a  9E  n 

V  •  B  -  -  r  •  =  0 

1  c  3t 


(47) 


_  -  1  a  9E  9B 

V.  x  E  —  r  x— —  - — 

1  c  9r  9t 


(48) 


_  -  1  a  9H  -  3D 

V,  X  H - r  X  — r—  =  1  +  — r — 

1  c  9t  j  9t 


(49) 


In  the  section  on  geometry  transformation  coefficients, 
such  that 


a  „  are  determined 


e.  =Z  a  . .  u.  (50) 

1  j  31  3 
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W  ‘Va'' 


Where  e  are  spherical  polar  unit  vectors  and  u.  are  prolate  spheroidal  unit 
i  3 

vectors.  In  particular 

r  =  #1 1?  +  <*21^  (51) 

Substitute  this  into  Eq.  (48)  and  (49)  to  obtain 

an  -»  _»  1  A  a  a  A  A  A 

+  v.  x  H  -  -  (a.A  +  ?  )  x  -r—  (H_  £  +H  t  +  H,4>)  (52) 

9t  J  1  c  11  21  9t  §  ?  9 

aR  -»  1  A  Aft  A  AA 

~=-V,xE+i  +„91?)xf  (E.?  +E  t  +  E,<t>)  (53) 

9r  1  c  11  21  or  €  ?  9 

Recall  that  the  prolate  spheroidal  unit  vectors  are  cyclic  in  order  ?  -£  -<j>, 
then  these  reduce  to 


_9D 

9t 


-7  +  V .  X  H  -  - 

J  1  c 


9H 


4> 


9H 


21  9t 


a 


11  9t 


a  A 

9  e 


9H 


9H, 


+  I  CL 


11  9  r  a2\  9t 


e 


(54) 


_9B 

9t 


•  V,  x  E  +  - 
1  c 


or, 


21  9t 


8E<,  a  %  a 


11  9t 


9  E. 


+  a 


K 


9E 


11  a21  9t 


«\  J 


(55) 


Assume  azimuthal  symmetry  and  write  the  remaining  terms  of  the  above 
equations  as  functions  of  ?  ,  £,<{>.  Also,  set 
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D  =  €E, 


B  =  mH 


(56) 


Then 


A 

€ 


A 

E  $ 

r 


E  $ 

99 


0 


-  c 


2e 

U  5 

T 


a  )(l  - 


T 
€  ) 


_9_ 

a? 


-  a2)(l  - 


§ 


\M 


2A 

u  K 

IT 


7“T 

a  )(? 


_ a_ 

“7) 35 


-  a2)(l  -  5 


+  U  <t> 


2A  ^ 


2  2 
a  Ml  -  ?  ) 


2  2  2 
?  -€  a 


9B 


Tf 


'A  A 

e 


9B 


or 


9 


11  9t 


A 

■? 


+ 


O 


21 


A 

9 


a  AAA 

—  (BSS  +B+). 


I  d 

i  rrr.  rr  "a? 


v(!2-eVki.I2) 


Mr2  2-  -2- 


(5  -  a  HI  -  C  )  E.  + 


A 

5 


(?2-  a2)(l  -  C2)  E 


A  v(?2  -  a2)(l  -  C2) 
9 


9 


2  2  2 
?  -5  a 


9 

a? 


-18- 


(57) 


+  UPJWMWmi  •W*-'**?**  MV# 


4a  -  »V  . \  j_ 

52  -  a2  V  '  85 


/  2  2  2 

v 

- T~  E? 

1  -  £ 


9E  A 

1  <(>  * 
+  —  O’  - —c 

c  21  8t 


3E.  A  /  3E  «E  \  A 

all  TT5  +  (  all  3t  '  a21  “37  )  * 


The  velocity  of  light,  u,  in  the  medium  is  determined  by 


Note  that  in  component  form,  Eq.  (57)  and  (58)  separate  into  two  uncoupled 

sets  of  equations  involving  E^  ,  E^,  ,  B^,  and  E^,  j^,  ,  and  B^  . 

We  assume  that  j  =  0.  Provided  that  the  second  set  of  fields  is  initially  zero, 
4> 

Maxwell's  equations  predict  they  will  remain  so.  Hence,  we  put 


WBS  =° 


The  field  equations  then  reduc  e  to 


if  + 


-?  a  )(1  -  ?  ) 


m2  -  a2)(l  -  C  2)  B, 


U  a21  n 


in  which  the  primed  fields  are  those  used  previously.  Maxwell's  curl  equa¬ 
tions  then  become 


9E_  j_  „  9B  2  9BX 

_ §_  °_F  +  _ 1  _ L 

9t  e  s  ~e  9?  ‘  c  9t 


(71) 


9  E 


t.  ° 


E. 


e  S 


i  „  9B,  2  9B , 

2  <j>  au  <(> 

—  -  u  a§  ‘  c  9t 


(72) 


9B 


<t 


„  2  2  9  E 

?  ~  a  £ 


9t 


2  2  2  0? 
5  -fa 


2  aE  2  9E 

1  -  f  f  a  1  - J  ? 

„2  .2  2  a?  '  c  2  2  2  9r 

S  -  £  a  S  -  £  a 


,  .2  2  9  E 

i  £  -  a _ £_ 

c.  ~2  2  2  dr 

C  -fa 


(73) 


b.  Boundary  Conditions  in  Retarded  Time 

In  the  limit  as  the  distance  between  two  points  approaches 
zero,  their  associated  retarded  times  become  equal.  Thus  intuitively  we 
expect  that  boundary  condition*  will  remain  the  same  in  retarded  time  as 
they  are  in  real  time.  This  is  uemonstrat-d  below  in  a  more  rigorous  man¬ 
ner. 

From  Maxwell's  differential  equations  (46)  through  (49),  equiva¬ 
lent  integral  equations  for  retarded  time  may  be  derived.  By  integrating  Eq. 
(46)  over  a  three-dimensional  space  volume  v,  we  obtain 


J  V.  Sdv^pdv  +  i  J 


a  9D 

r  *  ~Wr~  dv 


(74) 


Apply  Gauss's  divergence  theorem  to  the  left  side,  obtain 


/5  - 


D  dv 


(75) 


where  q  is  the  charge  in  v.  This  is  the  equivalent  of  Gauss's  law  in  retarded 
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time.  From  Eq.  (47),  a  similar  law  may  be  found  for  B 

/  B  '  dS  =^J rf*  Bdv 

s  v 

Thus  in  retarded  time  B  is  not  necessarily  solenoidal. 
From  Eq.  (48) 


(76) 


JvxE  •  ds  =  -J  de  +  ~  J$xE  ■  dB 
s  s  s 

Applying  Stoke 's  theorem  yields 


Je-  dt  -  + 


ds 


(77) 


(78) 


in  which  c  is  the  curve  that  bounds  the  surface  s.  This,  of  course,  is  the 
retarded  time  analog  of  Faraday's  law  of  induction.  Similarly  from  Eq.  (49), 
one  finds 


Jh  dl  =  J  j  ds  +  ip  j  D  •  ds  +  £  ~  Jr  x  H  •  ds 
c  s  s  s 


(79) 


which  is  a  generalization  of  Ampere's  law  in  retarded  time. 

Now  consider  the  boundary  conditions  in  retarded  time.  From 
Causs  law  Eq.  (75)  and  the  pillbox  shown  in  Figure  2,  one  finds 


/  5  •  d"s *  <Di  -  D2>  • ’-h*  kb ®i +  52>  ■ ?  ,80) 

si 

Here  it  is  assumed  that  D  may  be  discontinuous  on  the  boundary,  but  that  it 
is  continuous  in  regions  1  and  2  above  and  below  the  boundary.  The  lateral 
surface  of  the  pillbox  is  labeled  s^.  Take  the  limit  6 — ►(). 


***** 


where  X  is  the  surface  charge  density. 


m 


Region  1 


Region  2 


Figure  2 

Gaussian  Pillbox  for  the  Determination  of  Boundary 
Conditions  for  Normal  Field  Components 


Similarly,  from  Eq.  (47),  one  finds 


(B1  -  B2)  •  n  =  0 


From  Eq.  (76)  and  the  contour  shown  in  Figure  3,  it  follows  that 


(Ej  -  E2)  •  tAi 


1  A  “*■  A  A  1  0  ^ 


x  (E  +  E  )  •  n  x  tAj?6 

1  O 


The  contribution  to  the  line  integral  from  the  sides  6  has  been  omitted. 


Contour  Used  in  Determination  of  Boundary  Conditions 
for  Tangential  Field  Components 

Cancel  Aj 0  and  take  the  limit  6 — >0,  then 

(E1  -  E2)  •  t  =  0  (84) 

From  Eq.  (77),  a  similar  argument  gives 

-*■  -*■  A  -*  A  A 

(H  -  H2)  •  t  =‘k  •  n  x  t  (85) 

-*  A  A 

where  k  is  the  surface  current,  n  is  the  unit  normal  vector,  and  t  is  the  unit 

tangent  vector.  Equations  (81),  (82),  (84),  and  (85)  will  be  recognized  as  the 

boundary  conditions  in  real  time. 

For  problems  with  azimuthal  symmetry,  it  is  a  simple  matter 

to  show  in  real  time  that  B ,  vanishes  on  the  z  axis.  Let  us  look  at  the  de- 

<P 

rivation  in  retarded  time.  Integrate  Eq.  (79)  over  the  contour  shown  in 
Figure  4. 


I 

Figure  4 

Contour  for  Evaluating  B  on  the  z  Axis 


For  our  problem,  B  =  B ,  <j> .  Tne  azimuthal  symmetry  ensures  that  the 

<t> 

magnitude  of  B  is  constant  on  c.  Therefore, 

<2?r  r)  (86) 

r),  z 


J 


B  •  drj  *  B 


<t> 


Here  rj  is  a  coordinate  measured  perpendicular  to  the  z  axis.  Provided  rj 
is  small  enough  that  B  does  not  change  appreciably  over  the  surface  s,  we 
can,  for  points  on  the  z  axis  other  than  the  burst  point,  set 


A  “* 

r  x  B 


ds  =*■ 


9B 

9t 


<f> 


r?=0,  z 


/ 


sin0ds 


(87) 


Use  the  approximation 


sine  =* 


which  is  valid  for  rj  «z  -  a.  Then  Eq.  (87)  becomes 


(88) 
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cr(i,  i,  k- 1/2) 


1  r  .  .  >1  ff(  i,  ]»  * 

—  [e,  (i,  3*  k )  -  E  ( i,  j,  k-1  )J  - - 27 

At  L  s  ’ 


-,  3?(i.  3.  k-1/2)  u2 

|E?U  j,  k)+E?U  j,  k-1  )J  - - “  +  2A£ 

[b^U  j+1,  k-1  )  -  B0(i,  j,  k-1)  +  B^( i,  j,  k)  -  B0(i,  j-l,  k )] 

-A  tv*- j-  k>- v*-  j-  k'i)] 


[E?(i'  3»  k)  -  E?  (i,  j,  k-1  )]  = 


a(  i,  .1,  k-1/2  ) 
2e 


i  jca  3,  k-1/2)  u2 
[E^i,  j,  k)  +E?(i,  j,  k-1  )J - 7  2AC 

[b^U-1,  j,  k)  -  B^(i,  j,  k)  +  B^(i,  3>  k-1)  -  B^(i+1,  j,  k-1)] 

-37  [Vk  j-  k)‘ V‘- j'  k'u] 


/  „  2  2 
1  /  5  -  a 


j,  k)-B^U  3,  k-1)]  =2ML2  _c2a2j 

'  1, 3 

[E?(i,  j+l,  k-1)  -  E?(i,  3,  k-1)  +E?(i,  3,  k)  -  E?  ( 1,  3-1,  k )] 


These  equations  may  be  used  both  in  the  air  and  the  ground  meshes  provided 
appropriate  values  of  u,  cr,  and  e  are  inserted.  The  nonzero  field  components 
on  the  z  axis  can  be  obtained  from  equations  (31),  (32),  and  (36).  Equations 
(94),  (95),  and  (96)  are  then  solved  by  a  method  analogous  to  that  used  in  the 
real-time  code.  The  one  difference  in  the  methods  is  that  in  the  retarded 
time  code  the  fields  are  determined  in  the  air  mesh  along  a  line  of  constant 
j,  the  boundary  conditions  at  the  ground  are  applied,  then  the  fields  are  ob¬ 
tained  in  the  ground  mesh  along  this  same  line  in  the  direction  of  decreasing 
?  ;  whereas  in  the  real-time  code,  the  boundary  conditions  are  applied  last. 

To  facilitate  handling  the  boundary  conditions  at  the  ground, 
two  lines  of  points  are  defined  which  have  C  =0.  One  set  of  these  is  assumed 


to  be  in  air  and  the  other  to  be  in  the  ground.  The  mesh  near  the  boundary  is 
labeled  as  shown  in  Figure  1. 


The  boundary  conditions  at  this  interface  are  that  E^  and 

are  continuous.  This  implies  that  9E„  /9t  and  8B,/8t  are  continuous.  Thus 

S  9 

we  set 
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Mi  ,  j,  k-1/2)  3E  ( i  ,  j,  k-1/2)  [BE  (i  -1/2,  j,  k-1/2) 
K  a _ _  s  g _ s  a _ 

8t  SFt  2  8t 


9E  (i  +1/2,  j,  k-1/2)' 


8B.(i  ,  j,  k-1/2)  8B.(i  ,  j,  k-1/2)  f8B.(i  -1/2,  j,  k-1/2) 

9  a  9  g  1  9  a 


^(i  +1/2,  j,  k-1/2 


By  writing  Eq.  (97)  and  (98)  in  difference  form,  we  obtain 


,  r  -  cr(  i  +1/2,  j,  k-1/2) 

—  [E?(ia,  j,  k)-E?(ia,  j,  k-lj]=-— I - R - 

[*E  ( i  ,  j,  k  )  +  E  (i  +1,  j,  k  )  +  E  (i  +1,  j,  k-1  )  +  E  (i  ,  j,  k-1  )1 
^  ^  O  ^  O  ^  S  ^  S  -J 


J  <1  +1/2,  j,  k-1/2)  2 

-‘—Si - jsLVv^’-VV1’  j>k) 


+  Vv j-  k-‘>-  W1,  j>  k’1)J  -i^r[B*(ig+1' j>  k) 

-  B,(i  +1,  j,  k-1  )  +  B,(i  ,  j,  k)  -  B,(i  ,  j,  k-1  )] 

P  g  Pa  Pa  J 


a(i  -1/2,  j,  k-1/2) 


r  E  <  i  -1,  j,  k)  +  E  (i  ,  j,  k) 
L  5  3-  5  & 
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/ 


R 


2  2_ 

,  .  .  ,  .  v“l  1  __i _ I_L- —  |Ee(i  ,  j,  k) 

■EjY'1,  15  a 


r  -  a‘A5  /4 


-  Ep  <  i_>  j,  k-1)  +E?(ia-1.  j.  k>-E?(V1,  j'  k_1) 


4m  \  *  * 

S  a 


+  E.(i.,  J.  k)-Ec(i„,  j,  k-1 )  +  ( ig+l*  1.  k)-E|(ig+1’  J'  k_1)] 


V  g  ?  Y 


a  ,  LzMJ* —  [ 2E  ( i  ,  J,  k)  -  2E 

4C  at  2  2.-2..  L  5  a  * 

4CAT  g  -  a  A?  /4 


(i  ,  j,  k-1 )  +  E?  (ia-l,  ^  k) 


(i.-l,  j.  k-D  +  W1'  h  k>  '  W*’  k'°] 


S  a 


(100) 


Evaluate  3E  / 9T  at  (i  .  j»  k-1/ 2  ) 
?  a 


n  <y(ia,  j.  k-l/2  ) 

i_  [e5  (ia.  j.  k)  -  E?  (ia,  j.  k-1  )]  -  -  — ^ 


[E{<la'  j.  k)  +E?  (ia 


ie  <  V  )•  k'1/2) 

i  ,  j,  k-1 )  l  -  - - f  2A? 

a  J  n 


[B^(ia,  i+l.  k-1)  -  iyia,  j.  k-l)  +  B4(ia.  j.  k) 

H.k»]  -^K<V  i.k)-B,Ua.i.  k-1)] 


(101) 


Evaluate  9E?  /8t  at  (ig,  j.  k-l/2  ) 


l  r  1  <7<V  h  k'1/2) 

—  [E?(i  j,  k)  -Es(i  ,  j,  k-i)J  =  L 


2e 


[E?(i  ,  j,  k)+E?(ig,  j,  k-1)] 


j„(i  ,  j,  k-1/2)  2 


u 


g 


2AS 


[B*(i  .  j+l.  k-l>-B+(ig,  j.  k-l)  +  B^(ig,  j,  k)-B^ig,  j-1,  k)] 


"sv[VV j' k> '  W j-  k'n] 


g 


(102) 


Evaluate  Maxwell's  curl  equations  at  (i  +1,  j,  k-1/2) 

g 


,  ,  -i  Mi  ' H.  j,  k-1/2) 

i[E5(ig+l,  j,  k)-E5(iff+l,  j,  k-1)]  =  — S- 


2e 


r  1  k{ip+1'  j'  k_l/2) 

[E?(i  +1,  j,  k)+  Eg(i  +1.  j,  k-l)J  -  -  - 


+  2sr[B*(1g+1>  j+1>  k-1)-  W1' k-l)  + W1,  k) 


W1,  j_1’ k)]  -sbDW1, j- k)-  W1- j-  k_1)]  <103) 


1  r  1  a(ig+1'  j’  k"1/2) 

^[Es(ig+1,  j,  k)-E5(ig+l.  j,  k-l)]=-— S- 


2e 


r  j- (i  +l,  j,  k-1/2) 

[E?(i  +1,  j,  k)+E?(ig+l,  j.  k-1  )J  - 


X— 
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k)  - 


B  .( i  +1, 
<P  g 


k)  + 


B.(i  +1, 
<P  g 


j,  k-1) 


B  ,( i  +2,  j,  k-1  ) 
<P  g 


]-^[Vig+1‘ 3'  k)-B*V- k'n]  (104) 


2  2 

TrtW1'  j'  k>  '  B*<V1,  j'  k'U]  =  W  ?  2  . 


0 


E.U.+l,  1+1.  k-1)  -  E?(ig+1.  j,  k-1)  +E?(ig+1.  i,  k) 


r  g 


2 

,  .  1  1  -  ^  _  rE  (i  ,  j,  k)  -  E  ( i  +li  j>  k) 

E5(ig+1,  1-1.  k)J  -  2^^2-^rL  5l  g  J  5  8 


r  - 


2  2 

n  1  g  -  a 

+  E  (i  +1.  j.  k'1  >  ~  E?(ig+2’  ^  ^-1  J  cAt?2  _  a2^2 

3  g 


1  -  AC 


fE?(ig+l.  ).  k)  -  Ef  (ig+i.  j,  k-1)]  -  cAt  ?2  _  2 


e* 


rE5(lg+l.  1.  k)-E5(ig+l,  j,  k-1)] 


(105) 


For  given  j  and  k.  Equations  (99)  through  (105)  are  solved  simultaneously  for 
the  four  independent  field  quantities  at  the  air-ground  boundary  and  for  the 
three  field  quantities  at  the  point  (ig+l«  3>  k)- 
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4. 


GEOMETRY  FOR  THE  LOW-ALTITUDE  EMP  CODES 


By  referring  to  Figure  5,  define  prolate  spheroidal  coordinates  as 
S=oT<ri-r)'  Z  =l(r+r),  <t>  =  tan'1^  (106) 

&2L  1  £*  1  X 


Note  that  surfaces  of  constant  £  are  hyperboloids  of  revolution  and  those  of 
constant  ?  are  prolate  spheroids  (i.  e.,  prolate  ellipsoids  of  revolution). 

Also  observe  that  all  real  space  lies  within  the  range  of  coordinates 

?  a,  -1  <  £  <  1  (107) 


The  z  axis  is  separated  into  three  segments  determined  by 


x  =  y 


0  or 


/ 

£  =  l 

for  z  <  a 

(108) 

S  =  a 

for  -a  <  z  <  a 

(109) 

£  =-l 

\ 

for  z  <  -a 

(110) 

Burst 

Point 


Ground 

Surface 


Figure  5 

Geometry  for  the  Prolate  Spheroidal  Coordinate  System 


/ 

f 


The  x-y  plane  is  taken  as  the  ground-air  interface  and  according  to  Eq.  (106) 
is  the  degenerate  hyperboloid  determined  by  ?  =0.  The  burst  point  of  the 
weapon  is  on  the  z  axis  at  the  point  z  =  a  which  is  also  described  by  ?  =  a, 

5  =1.  The  wave  front  of  the  electromagnetic  pulse  generated  by  the  weapon 
spreads  out  in  air  from  the  burst  point  with  velocity  c  and  is  determined  by 
r  =  ct,  or  by 


?  -  a€  =  ct 


(111) 


Thus  for  fixed  t  the  wave  front  describes  a  straight  line  in  the  prolate 
spheroidal  coordinate  system. 

From  Eq.  (106)  we  find 


r:  =  ?  +  a? 
r  =  ?  -  a? 


(112) 

(113) 


Now 


2  2 
+  y  +  (z  +  a) 

(114) 

2  ,  .2 

(115) 

+  y  +  (z  -  a) 

From  these  it  follows  that  the  transformation  equations  from  cartesian  co¬ 
ordinates  to  prolate  spheroidal  coordinates  are 


€  \/x2  +  y2  +  (z  +  a)2  -  V^2  +  y2  +  (z  -  a)2  (116) 

?  v£2  +y2  +  (z  +a)2  +|  \/x2  +  y2  +  (z  -  a)2  (117) 


<)>  =  tan  1  — 
x 
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(118) 


where  principal  square  roots  are  intended.  The  inverse  transformations  are 


X  =v(?2  -  a2)(l  -  C2)  cos9 

(119) 

y  =  2  -  a2)(l  -  5  2)  sin 9 

(120) 

z  =?£ 

(121) 

Transformation  equations  from  a  spherical  polar  coordinate 
system  centered  at  the  burst  point  to  prolate  spheroidal  coordinates  may  be 
obtained  in  a  similar  manner,  we  find 

-  2(a/r  +  cosfl) 

*  1  +  \/l  +  (4a/r)(a/r  +  cos#) 

£  =  \  [1  +  VA  +  (4a/r)(a/r  +  cose)  J 

The  azimuthal  angle  is  the  same  in  both  coordinate  systems.  The  inverse 
transformations  are 


(122) 

(123) 


r  =  ?  -  a? 

-i A t  -  a\ 

9  =  008 

We  now  seek  the  scale  factors  h_  ,  h  ,  and  h,  such  that 

€  ?  9 

ds2  =  h.  2  d£  2  +  h  2  d?  2  +  h  2  d92 
5  C  9 


(124) 

(125) 


(126) 


f 

{ 

I 


where  ds  is  an  element  of  length.  To  find  these,  differentiate  Equations  (119) 
through  (121)  and  substitute  into 


ds 


2 


dx  +  dy  +  dz 


(127) 


Obtain 


.  2  t 
ds  =  — 


2  2 

■^di2+i 


1  -  € 


2  2 
5  -  a 


V  dt2+<e2 


a2)(l  -  €  2)  d4>2  (128) 


comparing  Equations  (126)  and  (128)  gives 


(129) 


(130) 


h 


<t> 


2  2 
-  a  )(1  -  5  ) 


(131) 


Because  Eq.  (128)  does  not  contain  cross  terms  in  d? ,  d? ,  and  d<£  it  follows 
that  we  are  dealing  with  an  orthogonal  coordinate  system.  Incidentally,  the 
coordinates  are  cyclic  in  the  order  §  -?  -<J>. 

By  using  Eqs.  (129)  through  (131)  one  easily  finds  expressions 
for  the  various  vector  operations  in  prolate  spheroidal  coordinates. 


Vf 


-  € 
-  S 


df_  A 

2  2  8C  § 
a 


_9f_  f 


+  \/(S 


2W1  e  2*  9f  A 

a  >d  -5 


(132) 
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Let  us  now  consider  the  vector  transformation  equations  be¬ 
tween  prolate  spheroidal  coordinates  and  other  orthogonal  systems.  For 
convenience  denote  the  base  vectors  in  the  prolate  spheroidal  system  by 
u.  and  those  of  the  other  system  by  e..  We  may  expand  e.  in  terms  of  the 

A 

U.. 

1 


A 

e 


j  3 


(136) 


But 


A  A 

e. 


.  •  e.  =6..  =  Z  (a..u.  )  •  (a*  .uj  =  2  a  a  (u  •  u  ) 
1  3  i]  k,i  ki  k  £j  £  kfi  ki  f]  k  £ 


-  Z  a,  .  a  .  6.  =  2  a,  .  a.  . 

kj  ki  jfj  k|  k  ki  kj 


(137) 


i.  e. , 


6.. 

ij 


(138) 


similarly  it  can  be  shown  that 


^ki^ii  ~  6ki 


(139) 


Equations  (138)  and  (139)  are  the  orthogonality  relations  for  the  transforma¬ 
tion  coefficients  Multiply  (136)  by  and  sum  over  i. 


Z/\  r  A  r  .  A  A 

>v£vaiiVfvruk 


(140) 


Equation  (139)  has  been  used.  Thus 
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A  A 

u.  =  2.  a.,  e. 
i  j  i]  3 


(141) 


Any  vector  may  be  expanded  in  terms  of  either  set  of  base  vectors. 


f  =2  f.e.  =  £  f.a..  u.  =  £  f' .  u. 
i  1  1  i»j  1  31  3  J  3  3 


(142) 


Here  f 1 .  are  the  components  of  f  in  prolate  spheroidal  coordinates  and  f. 

3  ^ 

are  those  in  the  system  with  base  vectors  e^.  Observe  that 


f\  =  £<*.. f. 
1  j  H  3 


(143) 


Multiply  this  by  a.,,  sum  over  i,  and  use  the  orthogonality  relation.  Obtain 

lik 


f.  =  Y  a..f’. 
1  f  31  3 


(144) 


An  easy  v  ay  to  determine  the  transformation  matrix  is  to  find  an  expression 
for  the  prolate  spheroidal  unit  vectors  in  terms  of  those  of  the  other  coordi¬ 
nate  system  and  compare  to  Eq.  (141).  Now, 


€  =  h  V€  , 


t  =  h?  VC  , 


$  =  hx  V4> 
9 


(145) 


we  wish  to  find  the  transformation  to  spherical  polar  coordinates  centered 
at  the  burst  point.  By  taking  the  gradients  of  Equations  (122)  and  (123)  and 
employing  Equations  (129),  (130),  and  (145)  we  find 


e  a  a  1  ^ 
5  =  -a—  r  - ~r  6 


(146) 


?  r  —ft 

hC  h5 


(147) 
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A 

Of  course  <t>  is  identical  to  the  corresponding  spherical  polar  unit  vector. 
By  comparing  Eq.  (141)  with  Eqs.  (146)  and  (147)  the  desired  transforma¬ 
tion  matrix  is  obtained 


(a.) 

ij 


The  inverse  matrix  is 


(a..) 

i] 


-1 


(148) 


(149) 


5.  ELECTROMAGNETIC  FIELDS  WITH  SPHERICAL  SYMMETRY 
a.  Finite  Difference  Code 

A  finite  difference  code  is  described  below  for  calculating  elec¬ 
tromagnetic  fields  for  problems  with  spherical  symmetry.  The  assumption 
is  made  that  the  conductivity  and  Compton  current  terms  are  functions  only 
of  the  radial  distance  between  the  source  and  field  points  and  time.  It  is 
well  known  that  source  distributions  of  this  symmetry  produce  no  radiating 
fields.  Here  this  occurs  because  the  symmetry  reduces  Maxwell's  equations 
to  a  particularly  simple  form  involving  the  radial  component  of  the  electric 
field  and  time.  Because  of  this  simplicity,  the  code  runs  very  rapidly, 
which  is  a  primary  advantage  of  the  method.  This  code  has  been  used  as  a 
test  case  for  the  two-dimensional  codes  described  above. 
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In  M.K.S.  units  and  spherical  polar  coordinates  Maxwell's 
curl  equations  are 


A 

r 


2  . 

r  sine 


r  9  9  1  «  raE 

-t—  (r  sin0E  )  -  — —  (r  E  )  + - — 

80  d>  0  r  sine 


r  sin©  I  94> 


IF  <r  sine  ) 


]4[k 


9E  n 

(rv  V 


9 

at 


A  A 
B  r  +  B J  +  B, 
r  9 


$] 


(150) 


A 

r 


2  . 

r  sine 


r-t- 

|_  do 


i  £  r9H 

i  +  - - 

r  sing  8tj> 


(r  sine  H  )  -  —  (r  H  ) 
<j>  94>  o 


dr 


(r  sine 


a  3H 

—  (r  H  )  -  — r 
9r  0  do 


A  A  A 


9  A  A  A 

+  w  <v  +  V  +  D+*> 


(151) 


Assume  the  constitutive  equations 


— ^ 

B  =»H, 


D  =  eE 


(152) 


with  n  and  e  constant.  Also,  because  of  the  assumed  azimuthal  symmetry 
no  physical  quantity  depends  on  4>.  With  this,  Eq.  (151)  and  (152)  yield 


A 

r 

2 

r  sin0 


[w  (rsinSE4.)]’FiiH5[^ 


(r  sin0E,  ) 
9 


A  A  A 

r  +  B„8  +  Bx  ^  ^ 

r  e  <j> 


(153) 


- T7"  *r  sin0  Ba  )  '  9C  T~  <r  sin0  Ba  > 

r2  sine  89  +  r  sln9  3r  * 


9  .  .  8  Br  ‘'r  A  je  A  A 

[l7  (r  V'  8Fj=Tr+r  9+T4, 


9  ,  A  A  A 

+  (E  r  +  E  +  E  <j> ) 
8t  r  6  <|> 


(154) 


where 


(155) 


Note  that  in  scalar  form,  these  reduce  to  two  sets  of  coupled  equations.  One 

set  involves  E,  ,  B  ,  B  .  and  i ,  .  The  other  relates  E  ,  E  ,  B,  ,  j  ,  and  j  . 

9  r  0  9  r  0  9  r  0 

We  assume  that  the  current  distribution  is  such  that  i ,  =  0.  Then  if  E  ,  B  , 

J<j)  4>  r 

and  B  are  initially  zero.  Maxwell's  equations  predict  that  they  will  remain 

6 

so.  Thus  we  set 


E  =  B  =  Ba  =  j ,  =  0 
4>  r  0  9 


(156) 


Equations  (153)  and  (154)  then  become 


8E 

r 

at 


(157) 


~T - I?  (rsin0IV 

r  sin0 


8E9 

at 


— ~t — —  ■%-  (r sin0  B, ) 
r  sin0  9r  0 


(158) 


,  9E 
1_ _ r 

r  9  0 


?^<rEe> 


(159) 


When  transformed  to  retarded  time.  Equations  (157)  through  (159)  provide 
the  basis  for  several  EMP  calculations  involving  azimuthal  symmetry.  If 
also  one  demands  that  no  quantity  depends  on  0  so  that  spherical  symmetry 
obtains,  it  follows  that 


(160) 


8E0 

at 


2 

—  ~(r  B.) 

r  9r  P 


(161) 


at 


~  (rE  ) 
r  9r  0 


(162) 


Observe  that  these  equations  are  decoupled  into  two  sets,  one  involving  E^, 

j  and  the  other  involving  E  ,  j  ,  and  B  Now,  because  of  the  azimuthal 
r  0  n  y 

symmetry  j  and  E  are  zero  on  the  z  axis.  The  spherical  symmetry  demands 
that  they  cannot  depend  on  0.  They  are  thus  zero  everywhere.  From  Eq. 

(162)  it  is  evident  that  B^  is  constant.  Thus,  if  B^  is  initially  zero  the 
above  reduces  to 
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X 


8E  j 
r  r 


(163) 


j  — +  <*E 
Jr  Jr  r 


(164) 


where  j  is  now  the  Compton  current  and  is  the  conduction  current.  Then 
Eq.  (163)  becomes 


8E  j  aE 
r  r  r 


(165) 


This  is  the  differential  equation  that  determines  the  radial  component  of  the 
electric  field  for  problems  with  spherical  symmetry. 

This  equation  may  be  put  in  difference  form  as  follows.  Let 


r  =  jAr, 


t  =  kAt 


(166) 


then  Eq.  (165)  may  be  expanded  about  the  point  (j,  k-1/2). 


[Er<j,  k>  -Er(j,  k-l>]  =  -i[jr(j.  k-1/2) 


+  u(j,  k-1/2)  Er<j,  k-1/2  )1 


jr(j,  k-1/2  )At 


E  ( j.  k )  = 
r 


Y  _  a(j,  k-1/2  )At 

\ _ 2e 

1  +a(  j,  k-1/2  )At 
2e 


(167) 


Er(j,  k-1  ) 


(168) 


Provided  the  current  and  conductivity  are  known  functions  of  time  and  the 
initial  field  0)  is  given,  the  field  E^r,  t)  for  a  given  distance  r  =  jAr 
and  for  times  t  =  kAt  >  0  may  be  determined  recursively  from  the  above 
equation.  Because  E(j,  k)  depends  only  on  the  field  E(  j,  k-1  )  are  not  at  the 
fields  at  other  values,  we  conclude  that  it  is  not  a  propagating  field. 


b.  Exact  Solution 

In  addition  to  the  finite  difference  code  for  calculating  fields 
with  spherical  symmetry,  a  computer  code  for  evaluating  the  exact  solution 
to  these  spherically  symmetric  fields  has  been  developed.  For  this  special 
physical  situtation  Maxwell's  equations  become 


V  x  H  =  j(r,  t)  +  o(r,  t)  E(r,  +  “at  =  ° 


(169) 


which  has  as  an  exact  solution 


t  r  t 

t)  =  -J  j(t')  exp 


CT(t")dt")  dt 


(170) 


The  code,  the  listing  for  which  is  in  Appendix  B,  is  for  numerically  eval¬ 
uating  the  above  integrals  for  arbitrary  currents  and  conductivities  of  the 
form  of  Eqs.  (1)  and  (2).  This  technique  might  be  of  value  for  obtaining 
the  fields  inside  of  300  meters  to  be  used  as  boundary  conditions  for  a 
finite -difference  code. 


See  Longmire,  C.  L. ,  Close-In  E.M.  Effects,  LAMS-3072,  E-3, 
May  1 96  4 . 


II.  DIFFUSION  EQUATION 


1.  PROGRAM  AND  CALCULATIONS 

A  code  for  the  CDC-6600  at  the  Air  Force  Weapons  Laboratory, 
Kirtland  Air  Force  Base,  has  been  written  and  used  to  calculate  EMP 
above  an  infinitely  conducting  ground  in  the  diffusion  approximation.  In 
the  general  case  of  Maxwell's  equations,  one  obtains 


V  B  =  -/iV  x  J  + 


32B 

8t2 


4- 


Vo  x  E  +  a 


9B 

at 


] 


(171) 


and  if  the  conductivity  o  is  a  slowly  varying  function  of  space  and  if 


mo 


2  — 

3B 

l 

a  b 

at 

*  2 

c 

at 

(172) 


as  may  be  the  case  in  a  highly  conducting  medium,  then  Eq.  (171)  reduces  to 

V2B  -no  4?-  =  -/uV  x  J  (173) 

at 

the  inhomogeneous  diffusion  equation.  In  the  particular  case  of  azimuthal 
symmetry  cylindrical  coordinates  are  useful  and  the  diffusion  equation 
becomes 


32b,  aBA  aj 

- 5-  -M  T  <«  -g£  --  -It  -jf 


8z 


(174) 


* 

which  can  be  solved  by  using  Green's  functions  . 
Graham,  W.  R. ,  Babb,  D.  D.,  RM-4905-PR, 
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A  complete  description  of 
January  1966. 
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the  code  and  its  listing  will  be  described  in  a  forthcoming  EMP  Theoretical 
Note.  Numerous  calculations  have  been  made  for  various  weapons  at  various 
distances.  Some  illustrative  cases  are  included  in  Appendix  III. 

2.  EARLY  TIME  FIELDS  AND  LIMITS  OF  VALIDITY 

When  the  conductivity  is  low,  at  early  times,  the  condition  Eq.  (172) 
is  not  satisfied  and  one  knows  that  the  diffusion  equation  approach  is  not 
valid.  Also  in  the  solution  by  Green's  functions  there  is  a  homogeneous  term 
(see  footnote  on  page  48)  that  involves  the  initial  field.  This  field  from 
the  earliest  times  (as  calculated  by  the  diffusion  equation)  is  clearly  incor¬ 
rect,  but  its  influence  was  felt  to  be  negligible  on  the  peak  values  of  the 
magnetic  field  because,  it  was  argued,  the  homogeneous  part  caused  by  the 
Compton  currents  was  so  extremely  large.  Numerous  runs  on  the  computer 
gave  evidence  that  this  was  not  always  the  case.  Hence,  analytical  examples 
that  showed  that  for  some  weapons  the  early  time  fields  were  unimportant  to 
the  peak  value  were  found,  but  for  other  cases  they  contributed  as  much  as 
half  the  peak  field.  That  is  not  to  say  that  even  in  the  most  unfavorable 
cases  were  the  peak  answers  in  error  by  50  percent,  but  rather  that  half 
the  contribution  was  in  error  by  some  undetermined  amount.  A  reasonable 
estimate  is  that  the  undetermined  amount  would  surely  be  less  than  the 
value  calculated.  The  net  consequence  of  this  result  is  that  in  some  cases 
one  should  attempt  to  do  better  on  the  early  time  fields. 

3.  CRITIQUE  AND  SUGGESTIONS  FOR  THE  FUTURE 

From  these  studies,  one  may  say  that  the  diffusion  equation  approach 
to  calculating  the  EMP  from  a  nuclear  burst  is  a  simple  and  rather  accurate 
method  provided  that  the  early  time  fields  are  small,  or  are  handled  in  a 
proper  manner,  and  that  the  ground  approximates  a  perfect  conductor.  The 
method  is  considerably  simpler  than  a  finite  difference  scheme  in  that  it 
takes  far  less  computer  time  and  is  not  plagued  by  such  questions  as  stability 


or  convergence.  One  problem  that  has  occurred  in  this  approach  is  the 

question  of  the  size  of  the  time  cuts.  The  time  dependence  of  the  current 

* 

and  conductivities  has  been  approximated  by  a  series  of  exponentials  of  the 
form 

0o{t  -  V 

J(r,t)  =  J(r,TQ)  e  u  u  ,  TQ  <  t  <  Tj  (175) 

0,<t  -  T  ) 

J(r,t)  =  J(r,  T  )  e  ,  T  <  t  <  (176) 

I 

I 

02(t  •  T2) 

J(r,  t)  =  J(r.  T2>  e  ,T2<t<T3  (177) 

.  ! 

.  i 


and  so  on  for  as  many  time  cuts  as  one  feels  are  needed  to  approximate  the 
given  current  curve;  and  similarly  for  the  conductivity.  It  was  found  empiri¬ 
cally  on  the  computer  that  the  values  of  the  calculated  fields  were  somewhat 
sensitive  to  the  size  of  the  time  cuts,  i.  e. ,  if  more  smaller  time  cuts, 

T  "  T  ,  T  -  T  ,  . . .,  etc.,  were  used  to  approximate  a  given  current  or 
1  U  «  1 

conductivity  curve  (say  20  instead  of  seven  or  eight),  the  values  of  the  fields 
calculated  changed  by  perhaps  10  percent  at  peak  value.  This  problem  was 
reduced  to  negligible  status  by  simply  increasing  the  number  of  time  cuts 
to  such  an  extent  that  the  field  values  were  not  changed  by  further  increases 
in  the  number  of  time  cuts.  Thus  it  was  thought  that  the  variations  observed 
were  due  to  increasing  accuracy  in  the  representation  of  the  current  and  con¬ 
ductivity  curves.  More  recently,  analytical  studies  of  the  question  have 
shown  that  it  is  not  a  matter  of  numerics,  but  rather  a  basic  feature  of  the 

Graham,  W.  R.,  Babb,  D.  D. ,  RM-4905-PR,  January  1966. 
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method  of  solution.  That  the  matter  is  not  a  very  serious  one  is  indicated  by 
the  convergence  of  the  numerical  values  with  the  increasing  number  of  time 

cuts  and  by  the  agreement  of  the  field  values  obtained  with  fields  calculated 

¥ 

by  other  approaches  . 

A  further  extension  of  this  general  method  to  nondiffusion  fields  and  to 
a  finitely  conducting  earth  has  been  found.  It  is  based  on  the  use  of  a  Green's 
function  not  for  the  diffusion  equation,  but  rather  for  the  complete  equation 


V2B  -  Me  ?-§  =  -M?  x  J  (178) 

8t  9t 

The  Green's  function  for  this  equation  is  well  known  if  a  is  a  constant 
The  present  effort  has  been  toward  extending  the  Green's  function  to  cases 
for  which  a  is  a  function  of  time.  This  has  been  accomplished  for  a  linear 
time  dependency.  In  the  future  it  is  hoped  that  this  work  can  be  extended 
to  more  useful  time  dependencies  such  as  an  exponential.  If  this  is  done, 
one  will  have  an  analytic  method  for  calculating  the  propagation  of  wave 
fields  in  conducting  media.  The  results  could  easily  be  checked  by  taking 
the  limiting  cases  of  very  large  a  (diffusion  equation)  and  a  =  0  (wave 
equation).  This  type  of  calculation  might  be  especially  useful  for  calcu¬ 
lation  of  close  in  fields,  or  for  a  layered  conducting  earth. 


As  for  instance  by  the  method  of  Graham,  The  Electromagnetic  Fields 
Produced  by  a  General  Current  Distribution  in  a  Conductive  Environ- 
ment  Under  Certain  Symmetry  Conditions,  Technical  Report  No. 
WL-TR-64- 1 53,  January,  1965. 


See  Morse,  P.  M. ,  Feshback,  H., 
Vol.  I. 


III.  GROUND  SOURCES 


1.  INTRODUCTION  TO  GROUND  SOURCES 

The  earth  under  a  nuclear  explosion  slightly  above  the  surface 
is  bombarded  by  neutrons,  causing  it  to  become  a  source  of  gamma  radi¬ 
ation.  This  radiation  in  turn  produces  Compton  current  and  ionizes  the 
air.  Because  these  effects --Compton  current  and  ionization  of  the  air--are 
sources  of  EMP  phenomena,  they  are  called  "ground  sources.  " 

The  intensity  of  the  radiation  sources  created  by  monoenergetic 
sources  of  neutrons  has  been  calculated  by  a  Monte  Carlo  computer  code 
run  at  the  Los  Alamos  Scientific  Laboratory.  The  results  of  this  run  have 
been  analyzed  to  separate  out  the  manner  in  which  gamma  sources  vary  with 
time,  range,  depth,  and  gamma  energy.  They  have  also  been  separated  into 
two  components  that  can  be  ascribed  to  different  kinds  of  neutron  interaction 
with  dirt,  namely  inelastic  collisions  and  capture  collisions.  The  results  of 
this  analysis  are  the  six  S  functions  described  in  Section  2-b  below  In  order 
to  make  the  ground  sources  continuous  with  respect  to  space  and  time,  four 
of  the  S  functions,  those  with  a  spatial  or  temporal  dependence,  have  been 
approximated  by  functions  that  are  continuous  with  respect  to  space  and  time. 

Gamma  radiation  from  the  earth  interacts  with  the  air  to  produce 
ground  sources.  The  ground  sources  caused  by  a  monoenergetic  gamma 
source  at  a  specified  point  in  the  ground  were  computed  by  the  Monte  Carlo 
computer  code  called  TIGRE.  The  results,  called  "TIGRE  quantities,"  are 
described  in  Section  2-c. 

We  have  worked  toward  combining  the  results  of  these  two  Monte  Carlo 
codes  in  such  a  way  that  they  are  usable  as  inputs  to  an  EMP  ground  burst  code. 

The  task  is  not  complete;  this  report  describes  the  current  state  of  the 
project.  Portions  of  the  task  reported  previously  are  not  redocumented  here. 
Final  results,  including  graphs  of  ground  sources,  will  be  reported  in  a 
memorandum  to  WLRPE. 
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2. 


CALCULATING  GROUND  SOURCES 


a.  Introduction 

INT  is  a  computer  code  that  calculates  ground  sources  pro¬ 
duced  at  some  point  (in  the  air)  and  time  (since  the  burst)  per  neutron  of 
some  energy.  The  point  is  described  by  its  range,  distance  from  the 
origin- -the  point  on  the  ground  directly  below  the  burst- -and  its  angle 
from  the  vertical.  Currents  are  described  by  the  components  J^,  which 

is  positive  for  a  current  going  away  from  the  origin,  and  J  ,  which  is 

9 

perpendicular  to  J  and  positive  for  a  downward  current.  Thus  an  upward 

O 

current  at  the  surface  would  be  represented  by  a  negative  J  (90  )  and  a 

9 

zero  J  (90°). 
r 

INT  always  calculates  ground  sources  for  a  standard  set  of 
11  neutron  energies  and  10  angles  and  always  computes  the  average  ioni¬ 
zation  rate,  Q,  for  each  neutron  energy.  For  an  input  range  and  time, 

INT's  standard  output  is  21  numbers--Q,  J  at  10  angles,  and  J  at  10 

r  9 

angles- -for  each  of  the  11  neutron  energies. 

3 

Q  has  the  units  ion  pairs /meter  •  second  ■  neutron;  the  Js  have 
the  units  amperes /meter2  •  neutron. 

The  calculation  of  ground  sources  can  be  thought  of  as  having 
three  parts.  The  first  is  the  calculation  of  gamma  sources.  The  second 
is  the  accounting  for  the  attenuation  of  these  sources  as  they  travel  through 
the  air.  The  third  is  a  mechanism  for  combining  the  first  two,  i.e.,  for 
obtaining  the  effect  at  some  point  in  the  air  caused  by  sources  distributed 
through  the  ground. 

The  sources  of  gammas  are  represented  by  the  function 

S(E.t,  E  ,  r,  z,  t),  which  has  been  expressed  in  the  form  S  =  S  .  •  S  .  •  S  . 

N  7  ei  ri  zx 

+  S  •  S  •  S  .  The  latter  six  functions  are  called  "S  functions.  "  They 
ec  rc  zc 

are  described  in  Section  2-b. 

The  effects  in  the  air  at  point  p  of  sources  in  the  ground  at 
point  g  vary  with  the  location  of  the  source  with  respect  to  the  burst,  its 
7  energy,  and  the  location  of  p  with  respect  to  the  source.  In  general,  an 
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effect  F  is  of  the  form  F(E__,  E  ,g,p,t)  =  T(E  ,  g,  p,t)*  S(E.T,E  ,g,t').  The 

IN  y  y  IN  y 

quantity  T  is  called  a  TIGRE  quantity  and  is  described  in  Section  2-c. 

An  effect  F  =  T  •  S  is  the  result  of  a  point  source  in  the  ground. 
It  is  necessary  to  take  into  account  all  of  the  ground  in  order  to  have  a  ground 
source  F(E^,  p,  t)  as  described  in  the  first  three  paragraphs  of  Section  2 -a. 

To  calculate 


the  ground  is  divided  into  range,  depth,  and  angle  bins  by  using  a  cylindrical 
coordinate  system  centered  on  the  burst  axis,  and  the  E^  spectrum  is  repre¬ 
sented  by  six  discrete  energies.  The  coding  which  accomplishes  the  indicated 
integration  is  described  in  Section  2-d. 

b.  Sources  of  Gammas 

Basic  Definitions:  The  distribution  of  the  secondary  source  of 
gammas  induced  in  the  ground  by  a  near -surface,  monoenergetic  neutron 
impulse  has  been  expressed  in  the  following  form: 

S(E„,E  ,  r,  z,  t)  =  S  . (E  ,  E  )  •  S  ,(E.T,r,t)  •  S  .(E..,z,t)  •  1010 
N  7  ei  N  7  ri  N  zi  N 


+  S 

ec 


(E  )  •  S 
7  rc 


%,r,t) 


S  (E  ,  z,  t) 
zc  N 


photon  energy  produced 
(volume  )(time)  (neutron) 


(179) 


Here  E^  is  the  energy  per  neutron  of  the  monoenergetic  neutron  source,  E^ 
is  the  energy  per  photon  of  the  secondary  gammas,  r,  z  are  cylindrical  co¬ 
ordinates  with  the  origin  at  the  surface  directly  below  the  point  of  burst,  and 
t  is  the  time  after  source  impulse. 
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Energies  are  in  MeV. 
The  units  of  the  six  functions  are: 


The  units  of  S  are  MeV/ meter 


3 


sec  •  neutron. 


S  . 
zi 


S 

zc 


S  . 
ri 


S 

rc 


dimensionless 


-1 

cm 


meters 


MeV  •  meters  •  shake 


MeV  •  meters  •  second 


The  "i"  terms  arise  from  inelastic  collisions;  the  "c"  terms 
arise  from  capture. 

The  four  functions  S  .,  S  .,  S  ,  and  S  have  been  curve  fit  as 

n  zi  rc  zc 

functions  of  the  spatial  and  time  variables.  E1VT  and  E  have  been  left  as  dis- 

N  7 

crete  variables.  EN  is  a  discrete  source  neutron  energy  designated  by  the 
"file  number"  or  "file.  " 

Four  of  the  functions  have  been  normalized.  Specifically: 


-oo 

-00 

IS  . 
E  ei 
7 

=  1, 

Is  =1, 

E  ec 

7 

i 

[  S  .dz  =  1, 

J  21 

0 

f  S  dZ  =  1 

J  zc 

0 

LASL  Bin  Definitions: 

The  six  functions,  before  any  curve 

fitting,  were  derived  from  the  output  of  a  Monte  Carlo  code  run  at  the  Los 
Alamos  Scientific  Laboratory.  The  range  of  each  variable  was  divided  into 
bins  as  follows: 
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Time  bins  in  microseconds 

0  -  3  -  10  -  30  -  300  -  3000  -  30000  -  300000 
Radial  bins  in  meters 

0  -  50  -  100  -  150  -  200  -  300  -  400  -  500 
Depth  bins  in  meters 

0  -  0.  1  -  0.3  -  0.  5  -  1. 

7  energy  bins  in  MeV 

0-1-2-3-5-8-12 

Files  2  through  12,  representing  11  monoenergetic  neutron 
sources,  were  used.  The  correspondence  between  file  number  and  neutron 
energy  in  MeV  is  as  follows: 


File 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

MeV 

14 

10 

8 

6 

4.  5 

3 

2 

1 

.  3 

.  1 

.  01 

All  files  assume  dry  ground.  The  pre-fit  values  of  the  six  functions  derived 
from  LASL  Monte  Carlo  data  are  listed  in  the  computer  output  hereafter  re¬ 
ferred  to  as  SOURCE. 


Functions  of  Gamma  Energy  (Se.  and  S^):  The  total  gamma 

energy,  E  ,  spectrum  has  been  divided  into  six  energy  bins  with  boundaries 
7 

at  0-1  -  2-  3-  5-  8-12  MeV.  S  and  S  .  give  the  fraction  of  the  total 

ec  ei 

energy  produced  in  each  energy  bin. 


S 

S 


ec  j 


.  ,  ,  ,  ,  ("inelastic  collisiorfl 

total  gamma  energy  produced  by  ( capture  j 


E 

7i+l 

E  dN(E  )  (181) 

7  7 

E 
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where  dN(E  )  =  number  of  j^aptur^e^}  Pro^ucec^  Photons  having  an  energy 

between  E  and  E  +  dE  .  The  values  of  S  .  and  S  used  by  program  INT 
777  ei  ec 

can  be  read  directly  from  SOURCE,  as  no  curve  fitting  was  done  for  these 
functions.  Note  that  S  is  the  same  for  all  files  since  it  is  independent  of 

G  C 

neutron  energy. 


S  :  S  has  been  curve  fit,  first  as  a  function  of  z,  then  as 
zc  zc 

a  function  of  time.  The  functional  form  assumed  is 


szc(z'  *> 


exp 


4/3^) 


-  exp 


z  -  2|32<t) 
431<t) 


21 


erf[ 


@3<t) 

+  exp 


04<« 


(182) 


where  z  >  0  above  ground  and  z  <  0  below  ground.  Here  (3  has  the  dimen - 

2  1 
sions  meters  ,  02  and  /3^  are  in  meters,  and  is  dimensionless. 


erf(x)  = 


VT-/e’y 

0 


dy 


(183) 


The  computer  code  that  fit  S  was  called  the  SZC  code.  For 

ZC 


each  time  bin  the  first  section  of  this  code  produced  betas.  The  bin  boundaries, 
z^,  were  greater  than  0  for  the  running  of  this  code.  A  corresponding  change 
was  made  in  the  assumed  form  of  S_  .  After  this  fit,  the  following  substitutions 


zc 


were  made:  =  1  /\Zj3^  and  0 =  l//3^.  Then  subroutine  TIMFIT  fit  each  set 

of  0  (i  fixed,  file  fixed,  time  bin  varying)  to  the  form  X  +  <*ln(t),  t  in  shakes, 


-57- 


with  the  following  time  segmentation: 


/3  (t)  =< 
Hri  ' 


X„,  +  ln(t)  3/u  <  t  <  3mseconds 
Mn  Mn 


X_  +  a „  ln(t)  3m  <  t  <  300mseconds 
Fn  Fn 


(184) 


for  n  =  1,  2,  3,  4. 

The  coefficients,  lambdas  and  alphas,  are  listed  in  the  document 
"Ground  Source  Terms  .  " 


S  S  .  has  been  curve  fit,  first  as  a  function  of  depth,  then 
zi  zi 


as  a  function  of  time.  The  functional  form  assumed  was 


S  .(z,  t) 
zi 


(185) 


X(t)  was  fit  to  the  form  X  +aln(t)  with  the  segmentation: 


\(t)  =  < 


XQi  +«1ln(t)  .06m  <  t  <  3/j  seconds 


XQ  +  aln(t)  3m  <  t<300M  seconds 


(186) 


The  dimension  of  lambda  are  cm.  The  argument  of  ln(t)  is  in  shakes,  and  z 
is  less  than  0  below  the  ground. 

The  values  of  the  lambdas  and  alphas  are  listed  in  the  document 
"Ground  Source  Terms.  " 


S^:  The  radial  dependence  of  the  capture  term  is  assumed  to 
have  the  form 


>:< 

Unpublished 


-58- 


Src<rit)  =  ~-^Z(t)  •  6(r)  +  [A(t)  +  B(t)  •  r]  exp£-r  •  C(t)]j  (187) 


The  following  assumes  that  S  is  slowly  varying  and  hence  can  be  taken 

zc 

outside  the  integral  of  the  total  source  function  and  fit  as  described  above. 

The  LASL  data  has  a  table  of  values,  S  rt(n,j),  for  each  file. 

rcO  J 

The  n  indicates  a  radial  bin,  the  j  a  time  bin.  The  problem  was  to  find  a 

function  S  (r,  t)  such  that 
rc 


r  T 

n+1  j+1 


SrcO<n-') 


'J  / 


S  (r,t)dt27rrdr 
rc 


(188) 


Ranee  Fit:  The  function 


1  f  Bl(j)  +  B2^  *  r  T  r  1 

Srcl(r-j)  =  2Vr  '  6(r)  +  - B^J) -  eJ“T  B^)J  '  (189) 


was  fit  to  the  SrcQ(n,j)s  by  a  range  fit  program  so  that 


a 

SrcO 


Srci(r»j)2,rrdr 


(190) 


for  each  j.  The  output  was  a  set  of  values  for  Bq,  B^,  and  B^.  A  few 
of  these  values  were  then  improved  by  hand  calculations. 

Time  Fit:  From  Eq.  (188)  and  (190)  it  can  be  seen  that  the 
desired  S^^r.t)  must  have  the  property 
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'  • 


(191) 


S^c2^r' t)  was  taken  to  be  S^^r,  j )/^\+j  “  ^j)*  To  simPlify  the  form  of  the  re¬ 
sult  the  following  substitutions  were  made: 


Z(j) 


Bo(j) 


B^j) 

A<j)  =  B3(j)<fj+1  "  Tj> 


(192) 


B(j)  =  Bg(j)(T.+1  -  T.)J  C(j)  =  B3(j)  (193) 


Thus 


Src2(r>  t)  =  ZjtT"  *  6(r)  +[A(^  +  B(^  *  r]  exp  [‘  CW  *  r]j  094) 

for  t  in  time  bin  j,  and  Eq.  (190)  is  satisfied. 

For  programming  purposes  S  was  separated  into  two  parts,  the 

rc 

first  corresponding  to  the  6(r)  term  and  the  second  corresponding  to  those 

terms  multiplied  by  the  exponential.  Neither  of  these  functions  of  time  was 

continuous.  Because  continuity  with  respect  to  time  is  highly  desirable, 

continuous  functions  were  constructed  in  the  following  manner:  Let  F(J)  be 

either  of  these  functions  for  some  fixed  r.  The  continuous  function  to  be 

constructed  shall  be  called  G(t).  Let  t  _  be  the  geometric  midpoint  of  time 

m  j  j 

bin  J  and  tj  the  upper  boundary  of  time  bin  J.  T^  =  0,  the  lower  boundary  oi 
the  first  time  bin.  G  is  to  be  constructed  for  time  bins  i  through  j. 
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3 


2  7  5F7  (A  +  B.  r)exp  [-  C  •  r] 


t  £  t  <  t  G(t)  =  F(i)("front  step") 
l-l  m,  i 


t  ,  j  ^  t  <  t„  G(t)  =  F(j)("back  step") 
3 


At  all  midpoints  t  T,  G(t  _)  =  F(J). 

m,  J  m,  J 


G(tj)  shall  now  be  defined  so  that  a  linear  interpolation  among 

these  points  and  the  G(t  T)s  will  satisfy 

m,  J 


i  t. 

r  J 

I  F(J)  =  / 

T  J 


G(t)dt 


This  requirement  is  slightly  less  strict  than  its  antecedent  Eq.  (191). 


front 

K-Ste^.{  continuous  center  — ••♦-back  stepVi 


f 

s 


To  satisfy  Eq.  (197)  we  define  G(t^),  see  illustration  above,  to  be  such  that 
(the  area  of  triangle  a)  =  (the  area  of  triangle  b).  That  is 


,  .  (tm,k+l  '  tk)  G(tm,k+1)  +  (tk  ~  tm,k)  GV 

k  "  (t„  -  t,.)  +  (t,_  -  t__ 


m,  k+1 


m,  k 


(198) 


for  all  k  from  i  to  j-1,  inclusive.  Equivalently, 


G(t,  ) 
k 


(t  ,  -  t.  )  F(k+1)  +  (t.  -  t  .  )  F(k) 

m,  k+1  k _ k  m,  k  , 


*m,k+l  Sn,  k 


(199) 


for  all  k  from  i  to  j-1,  inclusive. 
If  we  now  set 


T(0)  =  t 


i-1 


T(2)  =  t. 


T(N)  =  t.  (200) 
J 


T(l)  =  t  T(3)  =  t  .  , 

m,  l  m, x  1 


(201) 


then  we  may  write  the  interpolation  to  be  used  between  t  .  and  t  .  ,  as 

m,  x  m,  j-1 

Oft)  ■  o[t(M  -  1)]  +  (g[t(M)]  -  g[t(M  -  1)]]  T(M)  T(“m  (202) 

for  T(M  -  1)  <t  <  T(M).  ("continuous  center") 

To  compute  this  continuous  center  when  F  =  1/2 ir  Z  one  needs  a 
table  of  G[T(M)]  s.  These  were  computed  by  program  CFITC  and  are  found  in 
Program  INT  in  the  array  A Z. 

When  F  =  [l/(2?rr)]  (A  +  B  •  r)  exp£-  C  •  rj  the  g[t(M)J  s  depend 
continuously  on  r.  They  are  computed  for  program  INT  by  subroutine  SRC3. 
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A  listing  of  the  parameters  A,  B,  C,  and  Z  is  given  in  "Ground  Source  Terms" 
where  they  are  called  AA,  AB,  AC,  and  AZ,  respectively. 


S  The  functional  form  assumed  for  the  radial  distribution  of 
n 

the  inelastic  term  is 


S  .(r,  t) 
ri 


+  b2(t)  •  r  +  bg(t)  •  r2 


(203) 


Dimensions: 


MeV 

1  meter  shake 

.  MeV 

b2  2 

meter  shake 

.  MeV 

b3  ,  3 

meter  shake 


b„  meter 
4 


S  . 
ri 


MeV 

2 

meter  shake 


Most  range  fits  were  made  by  program  SRI,  third  edition.  For  a  description 

^  >!< 
of  the  range  fit,  see  "SRI  CODE  as  revised  by  Ralph  Day,  October  1967  .  " 

The  result  of  this  range  fit  is  a  function  continuous  with  respect  to  r  but  a 

step  function  with  respect  to  t. 


Unpublished 
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r 


Sri3(r,  J)  =  Sr.(r,  t)  =  [A(J)  +  B(J)  •  r  +  C(J)  •  r2J  exp  [-  r  •  R(J)]  (204) 


for  t  in  time  bin  J,  which  has  the  units  MeV/meter  shake.  A  continuous 
function  of  time  is  created  in  exactly  the  same  way  as  for  S  .  Let  F(J)  = 

O 

S>ri3^r'  J)  x  10  for  some  fixed  r.  The  description  of  the  creation  of  G(t) 
from  F(J)  on  pages  is  now  applicable  with  but  three  exceptions:  When 

F  =  Sri3  •  108  then  i  =  1,  j  =  4,  and  the  G  Jt(M)J  s  are  computed  by  subroutine 
SRI3.  The  units  of  G(t)  are  MeV/meter2  second.  A  listing  of  the  parameters 
A,  B,  C,  and  R  is  given  in  "Ground  Source  Terms"  where  they  are  called  BA, 
BB,  BC,  and  BR,  respectively. 

LASL  data  for  times  greater  than  0.  003  seconds  was  ignored  in 
computing  S  . 

c.  Gamma  Transport- -TIG RE  Quantities 

A  source  of  gammas  in  the  ground  has  the  functional  dependence 

S(E„T,  E  ,  r,  z,  t).  An  effect  F  caused  by  that  source,  be  it  a  J  ,  a  J_  ,  or  a  Q, 
N  y  r  u 

depends  on  the  strength  of  the  source  and  on  the  position  of  the  source  with 
respect  to  the  point  p  at  which  the  effect  occurs.  TIGRE  quantities  give  this 
latter  dependence. 

The  point  p  is  described  by  spherical  coordinates  with  the 
origin  at  the  surface  directly  below  the  burst  center.  The  point  g  in  the 
ground  at  which  S  is  evaluated  is  described  by  a  cylindrical  coordinate  system 
(r,  z,  a)  centered  on  the  burst  axis  with  z  =  0  being  the  surface  of  the  ground. 
Let  d1  be  the  angle  from  the  vertical  of  a  line  directed  from  g  to  p,  and  r'  be 
the  length  of  that  line.  Further,  let  z  be  the  vertical  coordinate  of  point  g. 
The  appropriate  TIGRE  quantity  is  Ta(E^,  z,  r',  d'),  a  =  q,  r,  or  9.  Thus 


Q(E  ,  E  ,  p,  g,  t)  =  T  (E  ,  z,  r',  d1)  •  S 
N  y  q  7 


J  (E  ,  E  ,  p,  g,  t)  =  T  (E  ,  z,  r',  d')  •  S 
r  N  7  r  y 
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(205) 


(206) 


-f $355.4^! 


VEN'Eyp'g-t)  =  Te(Eyz'r'',*')'  s 


(207) 


The  functional  dependence  of  S  is  not  specified  in  the  above  equa¬ 
tions.  An  effect  F  at  p  at  time  t  is  the  result  of  gammas  that  left  gr'/c  seconds 
earlier,  c  being  the  speed  of  light.  Thus  we  wish  to  calculate  S  at  time 
t'  =  t  -  r'/c.  Now  we  have 


Q(EN,E  ,p)glt)  =  Tq<E^,z,r',*')  •  S<EN,r,Eyz,t') 


J  (E  E  p,g,t)  =  T  (E  z,r ',*•)  •  S(E  ,r,E  z,t') 
r  N  7  r  7  N  y 


Wy**'*0  =  VEyz*r,>,)  ‘  s(EN’r’E7'M,) 


’(208) 

(209) 

(210) 


d.  Integration 

So  far,  a  ground  source  at  point  p,  caused  by  a  gamma  source  at 
point  g,  has  been  calculated.  To  find  a  ground  source  F  at  p  accounting  for  all 
the  ground,  one  must  integrate  through  the  ground.  To  this  end,  the  ground  is 
divided  into  bins  in  an  obvious  manner  by  the  (r,  z,a)  cylindrical  coordinate 
system.  F  is  computed  by  a  summation  over  these  bins  of  the  effect  caused 
by  a  source  at  the  center  of  each  bin  times  the  volume  of  that  bin. 

To  find  a  ground  source,  one  must  also  take  into  account  all 
gamma  energies.  Because  of  the  way  Sg.  and  Sgc  have  been  constructed,  a 
summation  suffices.  Thus 


F(EN’P,t)  -Jl  (T  •  S)  d(vol) 
g 

-Llll  T  •  S  •  (bin  volume)  (211) 

r  a  E  z 
7 
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For  each  of  the  six  E  ,  a  standard  set  of  seven  zs  is  used. 

y 

These  represent  midpoints  of  depth  bins  and  are 


E 

y 

z. 

Z2 

Z3 

Z4 

Z5 

Z6 

Z7 

0.5 

0.  01 

0.03 

0.05 

0.08 

0.1 

0. 15 

0.25 

1.5 

0.015 

0.05 

0.08 

0.  1 

0. 15 

0.2 

0.3 

2.5 

0.  015 

0.05 

0.08 

0.  13 

0.2 

0.3 

0.4 

4. 

0.  02 

0.06 

0.  1 

0.  15 

0.2 

0.3 

0.4 

6.  5 

0.  02 

0.06 

0.  1 

0.  15 

0.2 

0.35 

0.6 

10. 

0.  025 

0.08 

0.  14 

0.2 

0.3 

0.4 

0.6 

Thirty-six  bins  represent  the  variation  of  a  over  half  its  range. 
Because  of  the  symmetry  with  respect  to  the  plane  defined  by  the  burst  axis 
and  point  p,  the  remainder  of  the  variation  of  a  can  be  (and  is)  represented 
by  doubling  the  volume  of  the  bin. 

Range  bins  were  10  meters  wide  and  extended  to  500  meters  or 

to  the  greatest  distance  a  neutron  of  the  energy  under  consideration  could  have 

traveled  since  the  burst,  whichever  was  shorter.  A  term  for  r  =  0,  which 

arose  from  the  curve  fitting  of  S  ,  was  included. 

rc 

INT,  the  coding  described  in  Section  2-d,  has  the  number  and 
size  of  depth,  alpha,  and  range  bins  fixed  as  described  above.  This  has  the 
limitation  of  giving  rough  approximations  to  ground  sources  at  points  near 
the  ground.  When  the  point  at  which  ground  sources  are  being  computed  is 
close  to  the  ground,  then  all  of  the  ground  bins  appear  to  be  out  at  a  side  and 
none  of  them  are  directly  below.  For  "low"  ground  sources  a  different  grid 
should  be  used. 
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p 


grid  represents  ground  well  grid  represents  ground  poorly 

Figure  7 

A  Limitation  of  Rigid  Ground  Gridding 


The  TIGRE  quantities  giving  the  radial  and  theta  components  of 

current  assume  that  the  center  of  the  coordinate  system  is  directly  above  the 

point  gamma  source.  That  is,  they  assume  the  U,  V  coordinate  system  is 

Figure  8.  To  calculate  J  and  J.  in  the  Y,  Z  plane  the  following  transfor- 

r  u 

mations  must  be  used. 


Z 


Figure  8 

Coordinate  Systems  and  Variables 
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Assume  all  gamma  source  points  lie  in  the  X,  Y  plane.  Let  a 
typical  such  point  have  polar  coordinates  (p,  a).  The  range  for  a  is  [-ir/2,  tt/ 2J. 
The  X,  Y,  Z  coordinates  of  this  point  are  (p  cosa,  p  sina,  0). 

Because  of  azimuthal  symmetry,  it  is  sufficient  to  assume  that 
all  points  at  which  current  and  energy  information  is  required  lie  in  the  upper 
right-hand  quadrant  of  the  Y,  Z  plane.  Let  a  typical  such  point  have  polar 
coordinates  (r,  <f>).  The  X,  Y,  Z  coordinates  of  this  point  are  (0,  r  sin<j>t  r  cos <f>). 


Figure  9 
The  Y,  Z  Plane 


Figure  10 
The  U,  V  Plane 


Let  U,  V  be  the  axis  of  a  coordinate  system  having  its  origin 
at  (p,a),  its  V  axis  parallel  to  and  with  the  same  sense  as  the  Z  axis,  and  its 
U  axis  positively  directed  through  the  foot  of  the  normal  through  (r,  <t>)  to  the 
X,  Y  plane.  Let  (r',  $')  be  the  polar  coordinates  of  the  point  (r,  <t> )  relative 
to  the  U,  V  axes.  Then 

r'  = 
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Let  R  and  0be  unit  vectors  of  a  polar  coordinate  system  in  the 
Y,  Z  plane  as  shown  in  Figure  9.  Similarly  let  H'  and  0'  be  unit  vectors  of 
a  polar  coordinate  system  in  the  U,  V  plane. 

The  U,  V  components  of  a  unit  vector  directed  from  the  U,  V 
origin  toward  (r1,  are  (sin^1 ,  cos<£').  A  unit  vector  in  the  corresponding 
0  direction  has  the  U,  V  components  (cos<£',  - sin0T ) . 


J  =  J  ,R'  +  J.,0' 

r'  0' 


(213) 


be  a  current  vector  at  (r1,  <t>')  caused  by  a  gamma  source  at  ( p,.a ).  From  the 
above  considerations  we  see  that  the  current  in  terms  of  U,  V  components  is 


J  \  /sin cos^'N  /JrJ 


Jv/  \cos<t>'  "sin*7\J0. 


(214) 


Let  and  be  components  of  the  projection  of  J  onto  the  Y,  Z 


plane.  We  see  that 


JA  \  i 0 


o  1/ \ J 


(215) 


r  sin0  -  psina 

1  ~ 2  "  2  2 

’  p  -  2rp  sin0  sina  r  sin  <t> 


(216) 


j  is  cos|3  where  f3  is,  as  shown  in  Figure  8,  the  angle  between  U  and  Y. 


The  transformation  matrix  from  the  Y,  Z  system  to  the  R,  6 

system  is 


Thus 


(217) 


sin <f>  sin$'  +  cos <f>  cos <l>' 
cos $  sin -  sin$  cosf 


A  sin <f>  cos 0'  -  cos <f>  sintf' 


A  cos<A  cos<t>'  +  sin$  sin^' 
* 


(218) 


The  current  components  J^,  and  J^t  are  computed  by  taking  the 

product  of:  a  TIGRE  quantity,  T;  a  gamma  source  strength  per  unit  volume,  S; 

and  the  volume  of  the  bin  under  consideration.  See  Eq.  (211). 

INT  integrates  through  the  gamma  spectrum  and  through  the 

ground.  For  a  given  range  and  time  the  code  computes  for  each  of  11  neutron 

energies:  J  and  J  at  a  set  of  10  angles,  and  Q  averaged  over  these  angles, 
r  v 

The  angles  ai  e  the  midpoints  of  bins  9°  wide,  i.  e. ,  4.5°,  13.5°,  22.5°,  ..., 
85.5°.  It  has  been  run  for  ranges  2,  3,  5,  7,  10,  15,  20,  22,  24,  26,  30,  40, 
70,  100,  200,  400,  600,  800,  1000,  1400,  and  1500  meters  and  times  from 
0.3^seconds  to  2275.  ^seconds. 


3. 


USING  GROUND  SOURCES 


Ground  sources  are  desired  as  inputs  to  computer  codes  that  calculate 
the  EMP  effects  of  a  near-surface  nuclear  detonation.  More  specifically, 
they  are  required  for  the  code  G,  which  requires  ground  sources  at  a  sequence 
of  (range,  time)  points.  The  ionization  rate  is  required  to  be  a  single  number 
representing  Q  at  all  9  and  all  <t>,  for  the  specified  (range,  time)  point.  Cur¬ 
rent  density  is  required  to  be  broken  into  two  perpendicular  components, 

and  J  ,  of  the  spherical  coordinate  system.  J  is  required  to  be  analyzed  into 
9  r 

a  finite  series  of  ordinary  Legendre  polynomials;  J  is  required  to  be  analyzed 

U 

into  a  similar  series  of  associated  Legendre  polynomials.  In  both  series  only 
odd-numbered  polynomials  are  allowed.  Thus  G  requires  Jr  be  represented 
by  the  first  n  odd  coefficients  of  a  series  of  ordinary  Legendre  polynomials  and, 
that  J  be  represented  by  the  first  n  odd  coefficients  of  a  series  of  associated 

U 

Legendre  polynomials.  As  Q,  the  current  is  assumed  to  be  the  same  for  all  0. 

Ground  sources  are  computed  by  INT  on  a  "per -neutron"  basis  for  each 
of  several  discrete  neutron  energies.  Before  G  can  use  these  sources,  they 
must  be  converted  by  means  of  the  neutron  spectrum  to  numbers  representing 
all  neutron  energies. 

The  coding  that  integrates  over  neutron  energies  and  calculates  coef¬ 
ficients  of  Legendre  polynomial  expansions  of  functions  of  6,  is  called  GROUSE. 


INT  (RR,  TSEC) 


Computes  ground  sources  at  RR  meters  and  TSEC  seconds. 
Programmed  for  use  as  a  subroutine  in  a  restartable  program. 
Prints  inputs;  stores  results  on  TAPE  5. 


Ground  Zero  Component 


^✓ffi stance  frohv. 
ground  zero  to 
point  greater  than  largest' 
orange  for  which  TIGR^. 
pgSmi antities  knowi^X^ 


Call  TIMDEP,  i.  e. ,  set 
quantities  dependent  on 
time  only. 


Calculate  BZR. 


T  -  time  at  which  ground  zero 
sources  should  be  evaluated. 


^^Tless  thans\s 
least  time  at  which 
non-zero  sources 
appear  ? 


Do  125  J2-1.  NATB  air  6  bins  I 

Set  TIGRE  index.  It  depends 
on  air  6 . _ 


Calculate  SZC  at  each  of  the  42 
( y .  z)  comb inat ion s . 


CS  rJ  (J2)  =  BZR  £ 


SEC(J3)  ■  E 


JT  ^TIGRE  -  R 
K(z's)  '*T‘ 


TIGRE  | R)  (TIGRE  index. 


7  (J3).  z  (J3))  ■  SZC(zk(J3).  t,  L)  •  Azk<J3) 


Subroutine  SDATA 


- - n  Assume  R,  t,  L  have  been  set 

(y.  Enter  j  Return  y  sources  at  the  standard  y 
^  energiec  and  z's  via  array  SS 

Set  up  quantities  dependent  only  on 
time,  i.  e., 

Call  TIMDEP 
TLOG  = 


iSRI=SRI(L.  R.  t) 

311 - A — 


I  SRC=SRC(L,  R.t\ 

511 - i - 

fPQ  1  K=l.  61  energies  I 


i  SI =SRI-  SEl/E 

iSC=SRC- SEC/E. 


I  DO  1  M  =  l.  7  depths! 

/SRI  =  0 

]Ye8 

4 - 

'  r 

<^SRC  =  0  ?N-I 
ifes 


iS/SRC  =  0? 

Yno 

Cail  ErESZC,  i.  e.  ,  calculate  parfs~oT 
SZC  dependent  on  L  and  t  only 


SRI  =  0? 
|TNo 

Calculate  part  of 
and  t  only 


lependent  on 


1 - 1 

606  CONT 

SZI=SZI(L,  z(K,  M),  t)  I 

1 

SZC  =  SZC(L,  z(K,M)~0l 


SS(K,  M)  =  (SI-  SZI+SC.  SZC)-  DELPPTHlK,  M)-  AREA 

7  sources  volume  of  bin 

being  considered 
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APPENDIX  I.  EXAMPLE  PROBLEMS  FOR  THE 
FINITE -DIFFERENCE  CODES 


The  field  plots  shown  in  this  appendix  were  obtained  by  putting  the 
current  and  conductivity  pulses  of  equations  similar  to  Eq.  (1)  and  (2)  into 
the  real-time  prolate  spheroidal  code,  the  retarded-time  prolate  spheroidal 
code,  and  the  one -dimensional  code.  The  intent  is  to  compare  the  codes, 
thus  the  shape  of  the  current  and  conductivity  pulses  were  not  chosen  to  be 
representative  of  any  weapon.  The  burst  height  used  is  300  meters  and  the 
spatial  grid  size  is  five  meters.  The  plots  labeled  "3"  are  for  a  field  point 
35  meters  above  the  ground,  whereas  those  labeled  "4"  are  for  a  field  point 
15  meters  above  the  ground. 
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RETARDED  -  TIME  PROLATE  SPHEROIDAL  CODE 


6.06-07  8.0E-07  I.0E-06  I.2E-06  I.4E-06  I.6E-06 

LOCAL  TIME 


2.0E-07  4.0E-07  6.0E-07  8.0E-07  I.0E-O6  I.2E-06  I.4E-06  I.6E-06 

TIME 


ONE  — Dl  MEN  SI  ON  A 


REAL  TIME  PROLATE  SPHEROIDAL  CODE 


2.0E-07  4.0E-07  6.0E-07  80E-O7  I.0E-06  I.2E-06 

TIME 


RETARDE 


APPENDIX  II 


CODE  LISTING  FOR  INTEGRATION  OF  EXACT  SOLUTION 
TO  SPHERICALLY  SYMMETRIC  FIELDS 
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O  n  o  o  o  n  r>  r>  o  ooo  ooo  ooooooooo 


PROGRAM  FXSOL (OUTPUT ) 


TM  AX  -  UPPE"  TIMF  LIMIT  ON  INTEGRATION 
TR  -  P ISF  T IMr 

T®  -  LOCAL  TIMF  FOP  PEAK  OF  CURRFNT  PULSE 
TSLAG  -  TIME  LAG  OF  CONDUCTIVITY  PFAK  RFHIND  CURRFNT  PFA’K 
JP  -  PEAK  CUPP  ENT  AT  GI VFN  R 

TPS  -  LOCAL  TIMF  FOP  PFAK  OF  CONDUCTIVITY  PULSE 
Sp  -  PFAK  CONDUCTIVITY  AT  GIVEN  o 

PEAL  JZF°0, JP, INTI IOOO) 

DAT  A  R,TMAX, JZFRO,RGAMMA,SZERO,RRCTA, TP , TP, TSL AG ,C  ,  FPS/ 

1 240.,1.EC-6,6.1EB, 300. ,20., 200.  ,1 . 5E-7 , 1 .666F-7 , 1 . ?F-B ,3 . E8 ,  B.HSF 
2-12/ 

S FT  CONSTANTS 

RFPS  =-l  ./FPS 
TPS  =  TP  ♦  TSLAG 
TPSL  =  TPS  +  TP  ♦  4.E-R 
R2  =  R*R 
RC  =  R/C 

SOLVE  FOR  PEAK  CURRFNT  AND  CONDUCTIVITY  VALUES  AT  GIVEN  R 

JP  =  JZ  FRO/P  2*E  XP I -P /RGAMMA ) 

SP  *  SZER0/R2*EXO(-F/RBFTA) 

ALFA-S  AND  BFTA-S  FOR  SFCOND  AND  THIRD  FITS 

TEMP  =  AL0G(2.*JP) 

ALFAJ  *  TEMP/TP 
RFTAJ  =  TEMP/TP 
TFMP  »  AL0G(2.*SP» 

ALFAS  *  TEMP/TR 
BETAS  =  TEMP/TPS 

COEFFICIENTS  FOR  FIRST  FIT  TO  CURRENT 

TEMP  *  EXP( ?.*ALFAJ*(-TR)) 

TJ  »  l./U. ♦TEMPI 

SS  =  I AL  FAJ-ALF A J*TFMP  I/I1.+TEMPI**2 
41  =  SS/TJ 

R1  =  ALOGITJI-ITP-TRMAl 

COEFFICIENTS  FOR  FOURTH  r-IT  TO  CURRFNT 

TEMPI  =  FXP(RFTAJ*(TP+TR) ) 

TFMP  *  FXP(  2.*PFTAJ*TP) 

TF  =  TEMPI/I l.+TFMP) 

SF  =  TFMPI* ( BFTAJ-BE  TAJ* TEMP) /(I .♦TFMP) **? 

TO  =  TP*TP 

A4  *  SF*T()*TO/(  TF*TF  I 
R4  *  l./TF  -  A4/T0 

C  COEFFICIENTS  FOR  FIRST  FIT  TO  COMDUCTIVITY 


I 


* 


4 


TEMP  =  FXP I  ?«*ALFAS*I-TR)  ) 

TJ  *  1 ./« l.+TEMO) 

SS  =  (ALFAS-ALFAS*TFMP)/(l.+TEMP)Y*2 
41S  =  SS/TJ 

R1S  *  ALOM  TJ)-<  TPS-TR  )*A1S 
C 

C  COEFFICIENTS  FOR  FOURTH  FIT  TO  CONDUCTIVITY 

C 

TEMPI  =  EXP(RETAS*TPSI  I 

TFMP  =  EXPI  ?.*BFTAS*I  TPSL-TPSH 

TF  *  TEMPI/!  I. ♦TEMP) 

SF  =  TEMPl*(«eTAS-RFTAS*TEMP)/n.*TFMP)**? 

TO  =  TPSL 

A4S  =  SF*TO*TO/( TF*TF  I 
R4S  =  1./TF-A4S/T0 
DT  =  5./ I6.*C ) 

T  *  PC  ♦  PT 
N  =  ( TMAX- T ) /DT  +  1 
PRINT  35,R,C,FPS,TMAX,OT,JP,SP 
PRINT  3^ 

C 

C  FIND  VALUFS  Or  l  NT  F  Gl<  A  NO  S  FOR  TRAPEZOIDAL  INTEGRATION 

3 

no  io  i  =  l * n 

CALL  CURRFNTIT.TR, TP,  ALE  A J,RETAJ,RC*FR, Al,Rl  ,  A4.P4) 

CALL  CONOUCT(T,TP, TPS, TPSL .ALFA S,RFTAS,RC, SI G  ,A1S, BIS, A4S,R4S> 

CALL  TR AP  (  T, TR, TPS, TPSL, TMAX, ALF AS, BETA S,RC , S TGINT , A1S , Bl S , A4S ,  R4S 
1  ) 

INTUI  «  FR*FXPIRF°S*SIGINTI 
PRINT  40,T,FR,SIG, INTI  I ) ,SIGINT 
10  T  «  T+OT 
C 

C  no  INTFGP  AT  ION 

C 

SUM  *  0, 

NM1  «  N-l 
00  20  I  *  2, NM1 
20  SUM  »  SUM  ♦  INTI  I) 

AREA  »  0T/2.RI  INTIU*TNTIN)*2.*SUM> 

EFIELO  *  REPS*AREA 
PRINT  50,  A® EA, EFIELO 

r. 

C  OUTPUT  FORMAT  STATEMENTS 

C 

30  FORMAT! 1H0, l IX, *TI ME*  ,  1  5X , *CURPFN T* ,  1  OX , *CONDUC T I VI TY* ,9X , *  I  NT FGR A 
1N0  INTEGRAL  OF  SIGMA  BETWEEN  T  ANO  T  MAX*I 
35  EnPM  AT  1 1  HI ,  *P.  »*F7.2/*  C  «*E12.5/*  EPS  **F12.5/*  TMAX  =*E12.5/ 

1*  OT  -4E12.5/*  JP  =*E 1 2. 5/*  SP  **F1 2, 51 
40  FORMAT  I 5E20 • 5 ) 

50  FnRMAT! 1H0.YTHE  INTEGRAL  WITH  THE  AJWVE  SET  OF  TNTFGRANOS,  CALCULA 

1TED  BY  THE  TRAPEZOID  RULE  IS  *,F14.5/*  THE  F-F  l FLO  WITH  THIS  IWTEG 
2RATI0N  IS*fia.5) 

END 
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1 1  . . . 


SUBROUTINE  ''.URRFNT(T,TR  ,  TP,  ALFA  J,0F  TAJ,  RE  ,FR,A1  ,B1  ,A4,B4I 
C 

C  GIVFN  REAL  TIME  AT  PARTICULAR  R,  WILL  SUPPLY  CURRENT  FROM  ONE  OF 

C  FOUR  FITS 

C 

TU  =  T-Rf 

IF(TU-TP+T°  ) 60, 60, 5] 

51  IF(TU-TP)61,61,5? 

5?  IF(TU-TP-TP )6?,6?,6i 

60  FR  =  FXP ( A 1 *TU*B 1 ) 

RFTURN 

61  TT  =  TU-TP  +  TP 

FR  =  FXPULFAJ*TT|/M,+FXP(2.*ALFAJ*ITT-TRm 

RETUPN 

62  FR  =  EXP ( P  FTA  J*TU  )/< l.+FXPf 2.*BETA  J*{  TU-TPH) 

RFTUPN 

6?  FP  =  TU/I  A4+B4*T(J) 

RETURN  i 

FND 


i 


1 

j 


I 


< 


>r 


SUBROUTINE  CONDUCT  (  T ,  TP. ,  TPS ,  TPSL ,  At  FA  S,  BF  TA  S,  PC ,  S  TO ,  A  l  S  ,« IS  ,  A4S ,  04S ) 
IS » 

c 

r  GIVEN  REAL  TIME  AT  PAPTICULAR  R,  WILL  SUPPLY  CONDUCTIVITY  FRP*1  ONr 
C  OF  FIVF  FITS 

C 

TU  »  T-RC 

IFITU-TPS+TR 160,60,51 

51  IFITU-TPS161 ,61, 62 

52  IF  fTU-TPSL  162, 62, 65 

60  SIG  *  EXP ( A  1 S*TU+B IS ) 

RFTUftN 

61  TT  *  TU- TPS*TR 

SIG  «  FXP(ALFAS*TT>m.*EXP(2.*ALFAS*< TT-TRII) 

RETURN 

62  SIG  =  EXPIBETAS*TU)/(1.  +  EXP(2.*BETAS*(  TU-TPSM  ) 

RETURN 

65  IF(TU-1.E-5)141, 141,144 

141  SIG  *  TU/I  A4S+-B4S*TU) 

RETURN 

144  SIG  *  SIG/2. 

RETURN 

END 
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r>  r>  o 


SU8RUUT  I NE  TRAP  (  T «  TV  ,  T»S  ,  TPSL  ,  TMA  X,  ALFA S,  BETAS  ,KC  , S T5  I  NT ,  AIS  ,  Bl  S , 
1 A4S ,  R4S I 

TR APE70T  DAL  INTEGRATION  FOR  CONDUCTIVITY 

DIMENSION  SIG( 50 ) 

Tl  =  T 
TU  *  TMAX 
DT  =  ITlJ-TD/49. 

00  10  I  =  1,50 

CALL  CONOMC  T(  TL,  TR,  TPS  »  TPSL,  ALF  AS,  BFTAS.P.C,  SIGMA  » A1 S  »R1S » A4S, R*S  ) 
STG(I)  =  SIGMA 
10  TL  =  TL  ♦  OT 
SOM  =  0. 

DO  20  I  *  2,49 
20  SUM  =  SUM  f  SIG(  1 > 

SIGINT  =  0T/2.*(  SIG(l)*SIGf 501  +  2. *SUM) 

RETURN 

ENO 


//  EXEC  SFTHf.FQN 

/£ 
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EXAMPLES  OF  DIFFUSION  EQUATION  CALCULATIONS 
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